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EXPERIMENTAL PERFORMANCE EVALUATION OF A RADIAL-INFLOW 
TURBINE OVER A RANGE OF SPECIFIC SPEEDS 
by Milton G. Kofskey and Charles A. Wasserbauer 
Lewis Research Center 

SUMMARY 

An experimental investigation was made to determine the effect of specific speed on 
efficiency for a 4. 59-inch radial-inflow turbine. The range of specific speeds investi- 
gated (72 to 108) at equivalent design speed and pressure ratio was obtained by changing 
volume flow, based on rotor exit conditions. Chahges in volume flow were accomplished 
by the use of stators having throat areas nominally 50, 75, 100, and 125 percent of design. 
The turbine was operated with air as the working fluid. 

Maximum total and static efficiencies were obtained over the specific speed range 
of about 80 to 90. The peak total and static efficiencies were 0. 91 and 0. 87, respectively, 
for the 75-percent configuration. 

An understanding of the losses which contributed to the variation of turbine perform- 
ance with specific speed at design blade-jet speed ratio was made possible by an analysis 
which determined the magnitude of the various losses for each configuration. Stator 
loss was the predominant contributor to the decrease in efficiency as specific speed was 
reduced from a value of 86. Rotor incidence and viscous losses were the primary con- 
tributors to the decrease in performance when specific speed was increased above the 
value of 90. Stator exit static-pressure measurements showed that, at equivalent design 
speed and pressure ratio, rotor reaction increased as specific speed increased. 

Rotor exit total-pressure and flow angle surveys indicated that low losses were 
obtained near the hub region of the rotor for all configurations at equivalent design speed 
and pressure ratio. Comparatively high losses were obtained near the outer wall. These 
increased losses may have resulted from tip leakage effects with blade unloading, as well 
as from centrifugation of low- momentum fluid to this region. 



INTRODUCTION 


The current Brayton-cycle space-power technology program at the Lewis Research 
Center includes the experimental investigation of factors which influence the performance 
of small radial-inflow turbines. One such factor is the specific speed parameter, which 
relates the operating variables of turbine rotative speed, volume flow based on exit con- 
ditions, and ideal specific work to turbine geometry and aerodynamic performance. 

Reference 1 shows specific speed - efficiency correlations for a number of radial- 
inflow turbines of various sizes and for a wide range of inlet conditions. This reference 
shows that high efficiency is attainable for a specific speed range from 65 to 105, with a 
significant reduction in efficiency outside this range. However, turbine size, rotor tip 
clearance, and Reynolds number effects are present in the specific speed - efficiency 
correlations but are not examined separately. Therefore, the experimental investigation 
described herein was conducted to determine the specific speed effect on performance 
for a particular turbine size with rotor tip clearance and Reynolds number held constant. 

Two approaches were considered to achieve the range of specific speeds. One was 
to design and fabricate an optimized stator and rotor configuration for each specific 
speed point, and the other was to use several stators with one rotor configuration. The 
second approach was chosen because it would minimize time and cost of the program; 
however, less than optimum turbine configurations may have resulted, especially at the 
extremes of the specific speed range. 

The 4. 59-inch-tip-diameter radial-inflow turbine of reference 2 was chosen as the 
research turbine. The design specific speed for this turbine is 95. 6. Three additional 
stators having throat areas nominally 50, 75, and 125 percent of design were fabricated. 
The four configurations cover a specific speed range of 68 to 107 at equivalent design 
speed and pressure ratio. Each configuration was investigated over a range of turbine 
pressure ratios at equivalent design speed. 

This report presents the performance of the subject turbine for each configuration 
and shows the specific speed effect on turbine efficiency. Results are presented in terms 
of equivalent weight flow and efficiency at equivalent design speed over a range of pres- 
sure ratios. Internal flow characteristics are presented in terms of static pressure 
variation through the turbine and radial variation of exit flow angle and loss distribution 
at the rotor exit. 


SYMBOLS 

2 

g gravitational constant, 32. 174 ft/sec 

H' isentropic specific work based on total -pres sure ratio, ft-lb/lb 
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Ah specific work, Btu/lb 

J mechanical equivalent of heat, 778. 029 ft-lb/Btu 
N turbine speed, rpm 

N g specific speed, NQ^ 2 /(H’) 2 /^, ft^^/(min)(sec^ // ^) 
p pressure, psia 

Q volume flow (based on exit conditions), cu ft/sec 
Re Reynolds number, w/pr^ 
r radius, ft 

U blade velocity, ft/sec 

V absolute gas velocity, ft/sec 

Vj ‘ideal jet-speed corresponding to total- to static-pressure ratio across turbine, 

(2gJ Ah.^ 1 / 2 , ft/sec 

W relative gas velocity, ft/sec 

w weight flow, lb/sec 

cc absolute rotor exit gas flow angle measured from axial direction, deg 
y ratio of specific heats 

6 ratio of inlet total pressure to U. S. standard sea-level pressure, p'j/p* 

e function of y used in relating parameters to that using air inlet conditions at U. S. 

standard sea-level conditions, 0- 740 /y + l\ y/y-1 

y \ 2 / 

r\ & static efficiency (based on total- to static -pres sure ratio across turbine) 

^tot t 0 * 3 * efficiency (based on total- to total-pressure ratio across turbine) 

6 squared ratio of critical velocity at turbine inlet to critical velocity at U. S. standard 

L.I 2 

sea-level temperature, (V i/V* ) 

Cij X LI 

M gas viscosity, lb/(ft)(sec) 

v blade-jet speed ratio (based on rotor inlet tip speed), U^/Vj 
Subscripts: 

cr condition corresponding to Mach number of unity 
id ideal 

w outer wall 
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nail 


t tip 

1 station at turbine inlet 

2 station at stator exit 

3 station at turbine exit 
Superscripts: 

* absolute total state 

* U. S. standard sea-level conditions (temperature equal to 518. 67° R and pressure 

equal to 14. 696 psia) 


TURBINE DESCRIPTION 

The 4. 59-inch-tip-diameter radial-inflow turbine described in reference 2 was 
selected for this investigation. Air equivalent design values are as follows: 


Equivalent weight flow, we^^/5, lb/sec 0.616 

Equivalent specific work, Ah/0 cr , Btu/lb 11.9 

Equivalent speed, N^0 cr , rpm 29 550 

Equivalent total- to static-pressure ratio, P^/Pg 1. 540 

Total to total efficiency, 77 ^ 0. 880 

Total to static efficiency, « 0. 824 

0 

Blade-jet speed ratio, v 0. 697 

Specific speed, N , NQ^^/(H’ ft^^/(min)(sec^^) 95.6 

s 


The range of specific speeds at equivalent design speed and pressure ratio was ob- 
tained by changing volume flow by using stators with different throat areas. This was 
done by essentially changing the stator blade angle. Three additional stators having 
nominal throat areas of 50, 75, and 125 percent of design were used to obtain nominal 
specific speeds of 68, 83, and 107. 

Figure 1 shows the four stators and the rotor used in the investigation. The meas- 
ured stator throat areas were 49. 6, 75. 3, 96. 1, and 126. 1 percent of design. Here- 
after, each stator and rotor combination will be referred to as the 50- , 75-, 100- , and 
125-percent configuration. One stator blade of each configuration had an elongated 
leading edge to block the flow from entering the small area end of the inlet scroll. A 
description of the 100-percent configuration including velocity diagrams is given in 
reference 2. The 100- and 125-percent stators each have 14 blades, whereas the 50- 
and 75-percent stators have 18 blades each. In order to maintain acceptable stator-blade 
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75 Percent design 


50 Percent design 



100 Percent design 125 Percent design 

Figure 1. - Turbine rotor and four stators. 
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surface velocities, the 50-, 75-, and 125-percent stators have slightly different shapes 
than the 100-percent stator. 

It may be noted that, although the throat area of the 100-percent stator was 3. 9 per- 
cent smaller than design, results (as reported in ref. 2) showed that equivalent design 
weight flow was obtained at equivalent design speed and pressure ratio. Attainment of 
equivalent design weight flow results from the flow check procedure, in which the rotor 
throat area is increased by cutting back rotor trailing edges until equivalent design 
weight flow is obtained. 

The rotor has 11 blades and 11 splitter vanes. These splitter vanes are used over 
the initial third of the rotor, thereby increasing the solidity in this region. The resultant 
decrease in loading was required at the hub to prevent low blade pressure-surface gas 
velocities. 
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APPARATUS, INSTRUMENTATION, AND METHODS 


The test facility, instrumentation, and method of calculating performance parameters 
were the same as those described in reference 2, except that air was used as the working 
fluid. Figure 2 shows a cross-sectional sketch of the turbine test section and the instru- 
ment measuring stations. A varying area scroll was used to obtain uniform inlet con- 
ditions at the stator inlet. A center body was used at the rotor exit to obtain measurement 
of exit static pressure at the hub and at the outer wall. Radial surveys of total pressure, 
total temperature, and flow angle were made at the rotor exit. 

The 100-percent configuration was tested at nominal inlet conditions of 16. 0 pounds 
per square inch absolute and 540° R and resulted in a weight flow of 0. 657 pound per 
second at equivalent design speed and pressure ratio. A nominal Reynolds number of 
277 000 was calculated from this result; Reynolds number is defined herein as 
Re = w/jir^. In order to eliminate the effects of changes in Reynolds number on turbine 
efficiency, this parameter was held constant for all configurations at equivalent design 
speed and pressure ratio. Thus, the inlet total pressure was adjusted for the other 
configurations until a weight flow of approximately 0. 657 pound per second was obtained. 
Table I shows the values of inlet total pressure and temperature and the pressure ratio 



Figure 2. - Turbine test section and instrumentation. 


6 



TABLE I. - EXPERIMENTAL OPERATING CONDITIONS 


Configuration, 

Inlet total 

Inlet total 

Pressure- ratio 

percent design 

pressure, 

psia 

temperature, 

°R 

range 

125 

13.0 

536 

1.28 to 2.13 

100 

16.0 

540 

1.30 to 2.16 

75 

19.2 

540 

1.29 to 2.26 

50 

27.2 

542 

1.31 to 2.32 


range over which the turbine was inves- 
tigated for each configuration. 

The turbine was rated on the basis 
of both total and static efficiency. 
Turbine inlet and exit total pressures 
were calculated from weight flow, 
static pressure, total temperature, and 
flow angle. In the calculations of tur- 
bine inlet total pressure, the flow was 


assumed to be normal to the plane defined by station 1. The exit total temperature was 


determined from turbine power measurements. 


RESULTS AND DISCUSSION 

The results of this investigation are presented in two sections. The first section 
includes overall results in terms of equivalent weight flow and efficiency for a range of 
pressure ratios at equivalent design speed with cold air as the working fluid. The effect 
of specific speed on turbine efficiency is then shown. The second section discusses the 
internal flow characteristics of the turbine as determined from exit radial surveys of 
angle and total- and static-pressure measurements through the turbine at equivalent 
design speed and pressure ratio. 


Turbine Performance 

Weight flow. - Figure 3 shows the variation of equivalent weight flow we-^T^/S 
with inlet total- to exit static-pressure ratio at equivalent design speed. Equivalent 
weight flows of 0. 752, 0. 615, 0. 519, and 0. 367 pound per second were obtained for the 
12 5-, 100- , 75-, and 50- per cent configurations at the equivalent design pressure ratio 
of 1. 54. The variation of weight flow with increasing pressure ratio indicated that the 
flow was subsonic over the entire range of pressure ratios covered. The figure also 
shows that near choked flow conditions were obtained for the 50-percent configuration 
at the pressure ratio of 2. 32. The combination of near choked flow conditions obtained 
for the 50-percent configuration and the flattening of the weight-flow curves with decreas- 
ing stator-throat area indicates that the velocity level through the stator blade row was 
increasing with decreasing stator throat area. 

Figure 4 presents the variation of equivalent weight flow with stator throat area 
for equivalent design speed and pressure ratio. Equivalent weight flow is expressed as 
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Percent stator throat area 


Figure 4. - Variation of weight flow with stator 
throat area at equivalent design speed and pres- 
sure ratio. (Based on measured throat area of 
100-percent stator.) 


a percentage of the experimental equivalent weight 
flow obtained with the 100-percent configuration. 

The dashed line shown on the figure represents the 
case where equivalent weight flow is directly pro- 
portional to stator throat area. Comparison of the 
experimental curve with the ideal case shows that 
the weight flow increases at a lower rate than the 
rate of area increase. This indicates that the stator 
pressure ratio pg/p'j increased with increasing 
stator throat area and, therefore, rotor reaction 
increased. This change in rotor reaction resulted 
from the variation of stator to rotor throat area 
ratios of the four configurations. The change in 
rotor reaction among the four configurations is 
discussed further in the section Internal Flow 
Characteristics. 

Efficiency. - Figure 5 shows the variation of 
total and static efficiency with blade-jet speed ratio 
for each configuration. The highest efficiencies, 
at design blade- jet speed ratio, were obtained with 
the 75-percent configuration. Total and static 
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Blade - jet speed ratio, v 
(d) 50-percent configuration. 

Figure 5. - Variation of efficiency with blade - jet speed ratio at equivalent design speed. 

efficiencies were 0. 91 and 0. 87, respectively, for this configuration. These values are 
significantly higher than the total and static efficiencies of 0. 89 and 0. 83 obtained with 
the 100-percent configuration. At design blade-jet speed ratio, the lowest efficiencies, 
total and static, were obtained with the 125-percent configuration. These values were 
0. 85 and 0. 77 for the total and static efficiencies, respectively. 

The level of rotor exit velocity, as indicated by the difference between total and 
static efficiency, decreases with decreasing stator throat area. For example, at the 
design blade-jet speed ratio of 0. 697, approximately 8 points in efficiency are attributed 
to rotor kinetic energy for the 125-percent configuration, while only 3. 0 points in effi- 
ciency are attributed to rotor exit kinetic energy for the 50-percent configuration. This 
decrease in rotor exit velocity with decreasing stator throat area results from the change 
in the stator to rotor throat area ratio among the four configurations. The figure also 
indicates the variation of rotor exit kinetic energy with blade-jet speed ratio. Comparison 
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of figures 5(a) and (d) shows that there was a greater rate of change in exit kinetic energy 
with increasing blade- jet speed ratio (decreasing turbine pressure ratio Pj/p 3 ) for the 
12 5- per cent configuration than for the 50-percent configuration. This effect results from 
the larger variation of weight flow with pressure ratio for the 125-percent configuration 
than for the other configurations, as shown in figure 3 (p. 8). 

Figure 6 shows the variation of total and static efficiency with specific speed for all 
four configurations investigated. The dashed line representes the variation of efficiency 
with specific speed at the design blade-jet speed ratio of 0. 697. The upper plot in fig- 
ure 6 shows that the highest total efficiency value of 0. 91 was obtained at a specific speed 
of approximately 86. This efficiency value is 2. 0 points higher than the efficiency of 
0. 89, which was obtained at the design specific speed value of 95. 6 for the 100-percent or 
reference turbine configuration. It may be noted that the design blade-jet speed ratio 
curve (dashed line) passes through the peak efficiency point for all but the 50-percent 
configuration. The heavy curve shown in the figure represents the envelope of the effi- 
ciency curves for all configurations. This curve shows that maximum total efficiency 
is obtained in the specific speed range of about 80 to 90. 

The lower plot in figure 6 shows the variation of static efficiency with specific speed 
for the four configurations. The highest efficiency value of 0. 87 was also obtained at a 
specific speed of approximately 86. This value of efficiency is about 3. 0 points higher 
than that obtained for the 100-percent or reference turbine configuration at the design 
specific speed of 95. 6. The lowest peak static efficiency of 0. 77 was obtained at a 



Figure 6. - Variation of efficiency with specific speed at equivalent design speed. 
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specific speed of 111. It should be pointed out, however, that part of this decrease in 
static efficiency results from using the same rotor with each stator. Since the 125- 
percent configuration passes the largest volume flow of the four configurations, the rotor 
exit kinetic energy would be expected to be higher for this configuration. 

The variation of static efficiency with specific speed for design blade-jet speed ratio 
(dashed line) shows the same trend as the envelope curve represented by the heavy line. 
The highest efficiency of 0. 87 was obtained at a specific speed of 86, and the lowest effi- 
ciency of 0. 77 was obtained at a specific speed of 108. It may be noted that both total and 
static efficiencies obtained at design blade- jet speed ratio occur at or very close to the 
peak efficiency points for the 75- and 100-percent configurations and at lower values of 
efficiency for the other configurations. From these results, it appears that radial-inflow 
turbines should be designed for a specific speed range of about 80 to 90 for the attainment 
of high efficiency. 

Loss distribution. - In order to obtain an understanding of the losses which contri- 
buted to the variation of turbine performance with specific speed at design blade- jet speed 
ratio, an analysis was made to determine the magnitude of the various losses for each 
configuration. The method used involved the determination of velocity diagrams for each 
configuration from measured turbine work, weight flow, inlet conditions of pressure and 
temperature, speed, stator throat area, and results of rotor exit surveys of total pres- 
sure and flow angle. Design loss distribution between the stator and rotor was used to 
proportion the measured overall turbine loss for the 100-percent configuration. Stator 
losses for the other configurations were then assumed to vary in proportion to the average 
of inlet and outlet kinetic energy as determined from the velocity diagrams. 

Rotor incidence losses were determined through adjustment of the actual incidence 
angle, which resulted in an effective relative whirl velocity different from the velocity 
diagram value. The adjustment depends upon the blade speed, the number of blades, the 
rotor diameter, and the volume flow at the rotor inlet. The use of the effective relative 
whirl velocity is analogous to the use of the slip factor for centrifugal impellers. The 
remaining losses were attributed to rotor viscous losses. Figube 7 shows the results of 
these calculations. The various losses, expressed in terms of efficiency, are shown as 
a function of specific speed. The magnitude of the exit kinetic-energy loss is shown by 
the difference between total and static efficiency values obtained from figure 5 (p. 9) at 
design blade- jet speed ratio. 

Figure 7 shows that rotor incidence loss increases as specific speed increases 
above 90. This increase in rotor incidence loss results from an increase in the stator 
exit flow angle (as measured from tangential) and a decrease in the stator exit velocity 
with increasing specific speed. Rotor viscous losses also increase substantially with 
increasing specific speed. The increase in rotor loss results from the increased rela- 
tive velocity level through the rotor. Part of the increase in rotor viscous loss can be 
attributed to the manner in which the range of specific speeds was obtained. Figure 7 
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Figure 7. - Variation of turbine losses with specific speed at equivalent design speed 
and pressure ratio. 


shows that there is no significant change in stator viscous losses as the specific speed is 
increased above a value of 86. This would indicate that the combined losses resulting 
from the velocity level through the stators and the boundary-layer blockage did not change 
to any large degree. 

Decreasing specific speed below a value of 86 results in an increase in stator viscous 
losses. The figure shows that the losses increase from about 5. 0 points in terms of 
efficiency at a specific speed of 86 to about 11. 0 points at a specific speed of approxi- 
mately 72. These losses may be associated with the increased velocity level through 
the stator and the increased boundary -layer blockage due to a larger ratio of wetted area 
to flow area. 

Calculations also indicated that rotor incidence losses were insignificant below a 
specific speed of 90. Rotor viscous and exit kinetic-energy losses decreased with de- 
creasing specific speed, since specific speed is proportional to the square root of the 
exit velocity when rotative speed, rotor throat area, and pressure ratio are constant. 

Internal Flow Characteristics 

The determination of turbine internal-flow characteristics for each configuration 
was based on the measured static-pressure distribution through the turbine, together 
with the results of a radial survey of turbine exit total pressure and flow angle. 

Figure 8 shows the variation in stator exit static pressure with design stator throat 
area at equivalent design speed and pressure ratio for the four configurations investigated. 
Rotor reaction decreases and stator exit velocity increases with decreasing stator throat 
area, as was noted in the discussion of equivalent weight flow. 
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Figure 8. - Variation of stator exit static pressure with stator throat area at equiva- 
lent design speed and pressure ratio. 



Figure 9. - Variation of turbine exit flow angle with radius ratio 
at equivalent design speed and pressure ratio. 

The results of a radial survey of exit flow angle taken at equivalent design speed and 
pressure ratio are presented in figure 9 for the four configurations investigated. It may 
be noted that, as the stator throat area was reduced, the exit flow angle changed from 
predominately overturning (as denoted by negative angles) to underturning over the entire 
passage height. This trend in exit flow angle with configuration is to be expected since 
the rotor exit relative velocity decreases with decreasing stator throat area. 

The variation of exit flow angle and exit total and static pressure with radius ratio 
indicated that there was a nonuniform work distribution from hub to outer wall for all 
four configurations, with mi ni mum work occurring along the outer wall. This may be 
due to blade unloading which results from tip leakage and from centrifugation of low 
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momentum fluid to this region. 

Local values of total efficiency were 
calculated on the basis of the change in 
tangential momentum through the rotor and 
the radial distribution of total pressure at 
the rotor exit. These results are plotted 
in figure 10 in terms of turbine loss 
(1. 0 - as a function of radius ratio. 
The figure shows that the largest radial 
variation in loss or efficiency is obtained 
with the 75- and 100-percent configurations. 
However, the magnitudes of the losses for 
these two configurations are substantially 
lower than those for the other two config- 
urations. The curves show low losses 
along the hub region and comparatively 
high losses in the region near the outer 
wall for all configurations. 

Calculations were made from experimental results to determine the radius ratio at 
which the weight flow was divided into equal parts. A radius ratio of approximately 
0. 77 was calculated for all configurations. This coincides with the design mean stream 
line, as shown in figure 10. At the mean stream line the calculated local loss is approx- 
imately equal to the experimental value as obtained from overall performance for each 
configuration. 



Figure 10. - Variation of turbine loss with radius ratio at equiva- 
lent design speed and pressure ratio. 


SUMMARY OF RESULTS 

An experimental investigation was made to determine the specific speed effect on 
performance for a 4. 59-inch-tip-diameter radial-inflow turbine at equivalent design 
speed over a range of pressure ratios. Results are presented for operation at a nominally 
constant Reynolds number of 277 000 at equivalent design speed and pressure ratio. The 
effect of turbine size on performance was eliminated by use of the same rotor for each 
configuration. The range of specific speed values investigated at equivalent design speed 
and pressure ratio was obtained by changing volume flow through the use of stators having 
throat areas nominally 50, 75, 100, and 125 percent of design. From this investigation 
the following results were obtained: 

1. Comparison of actual equivalent weight flow with an equivalent weight flow, which 
is directly proportional to stator throat area, showed that there was a deficiency in 
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weight flow for the 125-percent configuration and a surplus of weight flow for the 50- and 
75-percent configurations at equivalent design speed and pressure ratio. This difference 
in weight flows was attributed to the corresponding changes in rotor reaction which result 
from the use of the same rotor for each configuration. 

2. Maximum total and static efficiencies were obtained in the specific speed range 
of about 80 to 90. In this range, peak total and static efficiencies of 0. 91 and 0. 87 were 
obtained with the 75-percent configuration at a specific speed of 86. The lowest peak 
value of efficiency was obtained with the 125-percent configuration. For this case, the 
total efficiency was 0. 85 at a specific speed of approximately 108. The corresponding 
static efficiency was 0. 77 at a specific speed of 111. 

3. An analysis of stator and rotor losses over the range of specific speeds investi- 
gated at equivalent design speed and pressure ratio showed the following: 

(a) Turbine losses were at a minimum for the specific speed range of about 80 to 90. 

(b) Stator viscous losses were the predominant factor in the decrease in total effi- 
ciency as the specific speed was decreased from a value of 86. 

(c) Rotor incidence and viscous losses were the predominant factors in the decrease 
in total efficiency as specific speed was increased from 90. 

4. Stator-exit static-pressure measurements obtained at equivalent-design speed and 
pressure ratio indicated that the highest rotor reaction was obtained for the 125-percent 
configuration. Rotor reaction decreased with decreasing stator throat area. 

5. Radial surveys of rotor exit total pressure and flow angle at equivalent design 
speed and pressure ratio indicated that minimum losses occurred near the hub region and 
the losses increased substantially near the outer wall. These increased losses may have 
resulted from tip leakage effects with blade unloading and from centrifugation of low- 
momentum fluid toward the outer wall. Comparison of losses between configurations 
showed that minimum losses were obtained from hub to outer wall for the 75-percent 
configuration. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, August 26, 1966, 

120-27-03-13-22. 
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mass of perturbating body, sun mass units 
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power, w 
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dynamic pressure, ~ p(V ! ) 2 , newtons/m 2 
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radius from origin to object, m 
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temperature, °K 

time, sec 

gravitational potential 

/ 2 

x,y,z accelerations due to gravity, m/sec 
absolute velocity, m/sec 
relative velocity, m/sec 
true anomaly, radians 

forces acting on object other than gravity, thrust, lift, drag, and 
perturbations due to perturbating bodies 

components of r, m 

angle between thrust and velocity vectors (sketch (a)), deg 
angle of rotation of thrust out of orbit plane (sketch (a)), deg 
power efficiency factor 
k 2 m r 

atmospheric density, kg/m 3 


oo argument of peri center, radians 

£ equatorial longitude of ascending node , radians 

Subscript: 

0 initial value 
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APPENDIX B 


VECTOR RESOLUTION 
Relative Velocity 

The relative velocity is defined as the velocity of the object with respect 
to the origin body. If the origin body is assumed to rotate about the z-axis, 
this velocity is given by 

V 1 = V - co X r (Bl) 

In x^y,z component form. 


7 x = V x + “V 

(B2a) 

v y - « 

(B2b) 

V z = V z 

(B2c) 


In the following sections, the atmosphere of the origin body is assumed to ro- 
tate as a solid body at the rate as. 


Thrust Resolution Along x,y,z Axes 


The thrust direction is specified with respect to the relative velocity vec- 
tor V 1 bv the angles a and 0 , as shown in sketch (a) (p. 4), Eor resolution 
of thrust vector into x,y,z compor^er^ts, it is convenient to define vectors^ A 
^nd P normal to and within the r,V 5 plane, respectively, such that V', A , and 
P form an orthogonal set. Thus, 


— > — > — > 

A = r X V = relative angular momentum per unit mass 

P = V' X A 

The thrust vector can then be resolved in the V* , A,P set as : 

P • V’ as FV' cos a 
F * A s FA sin a sin 0 
P * P = FP sin a cos 0 

Solving for ? yields 

— >■ F > > —> — > — > 

F = —x (V* cos a A X P + A sin a sin 0 P X V 1 + P sin a cos 0 P) 

-n£ 


(B3) 

(B4) 

(B5a) 

(B5b) 

(B5c) 

(B6) 
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or, in x,y,z component form, 

F x = X [v 1 cos o.(A y P z - AjPy) + A sin a sin p(P y V£ - P z V y ) + P sin a cos |B P x ] 

(B7a) 

F y = — [V cos a(A z P x - AjjPg) + A sin a sin p(P z V^ - P X V^.) + P sin a cos P P y ] 

p2 

(B7b) 

F z = p- [v 1 cos a(A x P y - AyP x ) + A sin a sin 3(P x V y - ByV^) + P sin a cos p P Z J 

(B7c) 


Aerodynamic Lift and Drag Resolution Along x,y,z Axes 

The drag vector D is alined with the relative velocity vector V' and is 
therefore given in x,y,z components as 


- V' V* V* 

D = -D — — - D -2- - D — 1 
V’ V* V' 


(B8) 


The lift vector 
fined orthogonal set 


may be resolved into components along the previously de- 
V' , A, and P by the following relations : 


L • V' = 0 (B9a) 

L • A = LA sin p (B9b) 

L • P = LP cos p (B9e) 

Solving for L yields 

L = -% (A sin p P X V* + P cos p P) (BIO) 

p 2 

or, in x,y,z component form, 

L x B i [a sin p ( P y V z - P z V y ) + P cos p P x ] (Blla) 

Pu 1— 

Ly = p [A sin p(P z V x - P X V’) + P cos p P y ] (Bllb) 

L z = p [A sin p(P x V y - P y Vi) + P cos p pj (Bile) 
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APPENDIX c 


TRANSFORMATION EQUATIONS PROM ORBIT ELEMENTS 
TO RECTANGULAR COORDINATES 



From spherical trigonometry used in reference to the celestial sphere shown 
in sketch (c), the following relations may be derived for the position coordi- 
nates : 


x = r(cos SI cos u - sin SI sin u cos i) 
y = r(sin SI cos u + cos SI sin u cos i) 
z » r(sin u sin i) 


where 


r =s 


and v can be obtained from 


u 


COS V SB 


1 + e cos v 
co + v 


cos E - e 
1 - e cos E 


and 


M = E - e sin E 


(Cla) 

(Clb) 

(Clc) 

(C2a) 

(C2b) 

(C2c) 

(C2d) 


The velocity components may be obtained by differentiating the position equations 

ii = = ~ v /^ and r = 

r 2 Vp 


using the two-body relations d = ir = ^ and r = t/- e sin v: 
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(C3a) 


where 


x = - (N cos i sin ft + Q cos ft) 

y = •*/— (N cos i cos ft - Q sin ft) 

▼ p 

z = (ET sin i) 


H * e cos co + cos u 
Q « e sin co + sin u 


(C3b) 

(C3c) 

(C4a) 

(C4b) 
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APPENDIX D 


RUNGE-KUTTA AND LOW- ORDER INTEGRATION 


SCHEMES WITH ERROR CONTROL 


The Runge-Kutta formula used is of fourth-order accuracy in step size h* 
It is of the form 


l 2 1 

s x 2 - X! *= -g (k! + 2k 2 + 2k 3 + k 4 ) 


X * a dependent variable 


X *= increment in the dependent variable 

J 1 

h 2 = increment in the independent variable t 
k-L = h 2 ^ 2 (t 1 ,X 1 ) 


k 2 - h 2 X 2 I X 1 + “o" 


k 3 = h 2 x 2 , X x + — 

k 4 = ^2^2 ^1 b- 2 y X x 4- k 3 ) 


A lower-order formula may be found by utilizing the three derivatives at 
t a tg, ti, and t 2 . If hj_ « tj_ - to and h 2 = t 2 - tj_^ the following Lagran- 
gian interpolation formula gives the derivative at any time tQ < t < t 2 ; 


. _ . (t - t-|)(t - to) 


X = Xq 


(t - t Q )(t - to) 


hi ( hi + h 2 ] 


. (t - t Q )(t - t x ) 
2 h 2 (hi + h 2 ) 


Integration of this equation from tj to t 2 yields 


,f _ i M 
A ' 6 W ) 


ho\ / -h. 


Xq + rp (h 2 + 3h]_)X]_ + (2h 2 + 


no 

1 + t -=■ 

hi 


hi X 2 


I 


The difference in the increments over the interval hg between the Runge-Kutta 
scheme and the_low-order scheme may be divided by a nominal value of the depend- 
ent variable X to obtain the relative error 52* Thus, 


S 2 - 


i - *] 


X 


(D4) 


The error is expected to vary as approximately the fifth power of h, which 
leads to 


5 » Ah 5 

(where A is a suitable coefficient) or in the logarithmic form 

log 5 * A 1 +5 log h 


(D5a) 


(D5b) 


where 


A' « log A (D6a) 

Let it be assumed that A 1 will vary linearly with t, the variable of integra- 
tion* Then A 1 at a time corresponding to tg can be found from A' at two 
previous points t^ and tg as 


AX - A' 

A-h ~ A-2 + 1 ( t 3 " *2^ (D6b) 

z 2 ~ Z 1 

and if hg = (tg - tg) and hg = (tg - tj_)/ 

+ <-*2 - A i > (“<=) 

Ct 


and on this basis S3 would be predicted to be 

log S3 A3 + 5 log h-2* 

It is desired that S3 should approximate the reference error > 

log I13 — *3 (log B* - A^) 


(D7) 

therefore, 

(JOB) 


Each dependent variable has an associated relative error and would lead to com- 
putation of a different step size for each variable^ however y the maximum rela- 
tive error of all variables may be selected for 8. Obviously^ inaccurate pre- 
dictions of step size can occur when the maximum relative error shifts from one 
variable to another or when any sudden change occurs. When a step size produces 
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an excessively large error (5 > 5 i imi t) r a reduced step size must be used. It 
may be obtained from the reference error ^ as 


h 3 



(log ^ - 



(D9) 


Starting the integration . - The Runge-Kutta scheme is simple to start, since 
integration from X^ to X^^. requires no knowledge of X less than x n . 

Since the error control coefficient A has no value at t = 0, a prediction of 
the second step size is difficult. To overcome this difficulty, two equal size 
first steps. may be made before checking the error. The A for the first step 
may be arbitrarily set equal to the A for the second step so that h 3 may be 
predicted. The low-order integration scheme equation in this case becomes, with 
h 2 = h]_, 



5 


(Xb + 4*l + ± z ) 


(DIO) 


Failures . - Should two consecutive predictions of the same step fail to 
produce an error 5 less than S i imi t.? a return to the starting procedure will 
be made with a third prediction on step size, which is no larger than one-half 
of the second estimate. The step-size control described here will operate stably 
with nearly constant error per step only for a well-behaved function. For most 
problems it will repeat a step occasionally to reduce a large error, and on sharp 
corners it will restart. This action is not regarded as objectionable. The ob- 
jective is to attain a desired level of accuracy with a minimum total number of 
steps. 
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APPENDIX E 




GLOSSARY OF VARIABLES 

VARIABLE 

COMMON 

LOCATION 

DEFINITION 

A (700) 

C(ll) 

ARRAY CONTAINING THE INITIAL OATA AND THE PROGRAM 
CONTROL VARIABLES 

A1 

B ( 10 ) 

ERROR CONTROL PARAMETER DEFINED BY EQ. (D6A) AT T(l) 

A 2 

BUD 

ERROR CONTROL PARAMETER DEFINED BY EQ* ( D6A ) AT T(2) 

AC0EE1 

B ( 12) 

INTERPOLATION POLYNOMIAL COEFFICIENT FOR VARIABLE 
STEP SIZE* EQ. (D3) 

AC0EF2 

B ( 13) 

INTERPOLATION POLYNOMIAL COEFFICIENT FOR VARIABLE 
STEP SIZE, EQ. ( D3) 

AC0EF3 

B ( 14) 

INTERPOLATION POLYNOMIAL COEFFICIENT FOR VARIABLE 
STEP SIZE, EQ. (D3) 

AEXIT1 (10) 

A( 103) 

ENGINE EXIT AREAS FOR AT MOST 10 STAGES, M**2 

AEXIT 

B ( 3 ) 

AEXIT1 ( NSTAGE ) 

AK (3) 

A151) 

RUNGE KUTTA COEFF IC IENT S, SET IN STDATA 

ALPHA 

A ( 49 ) 

ANGLE BETWEEN VELOCITY AND THRUST VECTORS, SEE SKETCH 
(A) 

ALT 

A ( 4) 

VEHICLE ALTITUDE ABOVE EARTH, M 

AM 

B ( 90 ) 

TOTAL VEHICLE ANGULAR MOMENTUM PER UNIT MASS, M**2/ 
SEC 

AMASS (30) 

A ( 347 ) 

PERMANENT LIST OF BODY MASSES IN ORDER OF PNAME LIST, 
SET IN STDATA, MASSES FROM ELIPS DATA BEGIN AT AMASS 
(21), SUN MASS UNITS 

AMC (3) 

B ( 87 ) 

X , Y , Z COMPONENTS OF ANGULAR MOMENTUM PER UNIT MASS, 
**2 ) /SEC 

AMSQRD 

B ( 91 ) 

SQUARE OF TOTAL ANGULAR MOMENTUM PER UNIT MASS,M**4/ 
SEC**2 

AREA! (10) 

A ( 1 13 ) 

AERODYNAMIC REFERENCE AREAS FOR AT MOST 10 STAGES, 
M**2 

AREA 

B(6) 

AREA1 (NSTAGE ) 

ASYMPT 

A( 7 ) 

SEE TABLE II 

ATMN 

A ( 2 1 ) 

SEE TABLE II 

AU 

A ( 29 ) 

ASTRONOMICAL UNIT,M 

AW (4) 

A( 55 ) 

RUNGE KUTTA COEFF IC IENTS, SET IN STDATA 

AZI 

A { 35 ) 

INITIAL AZIMUTH ANGLE, USED WHEN IMODE - 4, SEE SKETCH 
(B), DEGREES 
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8 (800) 


C(llll) 

ARRAY CONTAINING INTERNAL PARAMETERS NOT UNDER USER 
CONTROL 

BETA 


AI50) 

ANGLE BETWEEN VELOCITY-THRUST PLANE AND ORBIT PLANE, 
SEE SKETCH(A) 

BMASS ( 

8) 

BI137) 

BODY MASSES SELECTED FROM AMASS LIST IN SEQUENCE COR- 
RESPONDING TO BNAME LIST 

BNAME ( 

8) 

B ( 122 ) 

ORDERED LIST OF BCD BODY NAMES 

BODYCD 

(10) 

A ( 143 ) 

ORIGINAL UNORDERED LIST OF BCD BODY NAMES READ IN AT 
INPUT 

BODYL ( 

10) 

B ( 153 ) 

AUXILIARY ORDERED LIST OF BCD BODY NAMES 

GO 


A ( 1 65 ) 

TOTAL DRAG COEFFICIENT 

CD I 


A( 163 ) 

INDUCED DRAG COEFFICIENT 

CHAMP 


B ( 25 ) 

SMALLEST CRITICAL RADIUS (RBCRIT(J)) WITHIN WHICH 
OBJECT LIES 

CINCL 


B ( 55 ) 

COSINE OF INCLINATION 

CIRCUM 


B ( 82 ) 

CIRCUMFERENTIAL COMPONENT OF TOTAL PERTURBATIVE AC- 
CELERATION, M/ SEC* *2 

CL 


A ( 164 ) 

LIFT COEFFICIENT 

CLEAR 


C ( 3 ) 

SEE TABLE II 

COEFN ( 

192) 

A ( 407 ) 

STORAGE ARRAY FOR COEFFICIENTS USED TO COMPUTE ALPHA 
CL , CD I , CD OR OTHER PARAMETERS 

COMPA ( 

3) 

B ( 63 ) 

COMPONENTS OF TOTAL PERTURBATIVE ACCELERATION ALONG i 
, Y » Z AXES. 

CONSU 


A ( 31 ) 

SEE TABLE II 

CQNSTU 


A ( 32 ) 

SEE TABLE II 

COSALF 


B ( 48 ) 

COSINE OF ALPHA 

COSBET 


B ( 49 ) 

COSINE OF BETA 

COSTRU 


B ( 53 ) 

COSINE OF TRU 

COSV 


B ( 57 ) 

COSINE OF THE ARGUMENT OF LATITUDE 

D (1100) 

C ( 2 1 1 1 ) 

ARRAY WHERE SAVED DATA IS STORED FOR LATER USE. 
ARRAYS A, XPRIM,AND XPRIMB MAY BE SAVED. 

DELMAX 


A ( 19 ) 

SEE TABLE II 

DEL 


A (43) 

OUTPUT CONTROL PARAMETER USED IN STEP 

DELTl (10) 

A ( 1 33 ) 

INITIAL STEP SIZES FOR AT MOST 10 STAGES.SEC 

DELT 


am 

DELTl ( NSTAGE ) 

DNS I TY 


B ( 29 ) 

ATMOSPHERIC DENSITY, KG/M**3 

DONE 


8(39) 

CONTROL PARAMETER FROM STEP WHICH INFORMS NBODY TO 
STOP INTEGRATING 
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DRAG (3) 


DTOFFJ 
E 2 

EFMRS IT) 
ELEV 

EL IPS (12, 

EMONE 

END 

EPAR 

EREF 

ERLIMT 

ERLOG 

ETOL 

EXITA 

EXMODE 

FILE 

FLOW! (10) 
FLOW 

FORCE (3) 

GASFAC 

GEOH 

GK2M 

GKM 

H2 

1 BODY (8) 
ICC (10) 
IDENT (10) 

IMOOE 
IND (3) 
INDERR 


B ( 69 ) X,Y,2 COMPONENTS OF THE DRAG ACCELERAT ION, M/SEC**2 


A ( 23 ) 

B ( 18 ) 

B( 130 ) 

A ( 36 ) 

10) A ( 167 ) 

B ( 28 ) 
A( 5 ) 

B ( 26 ) 

A ( 13 ) 

A ( 14) 

B ( 17) 

A ( 30 ) 

B ( 392 ) 
B ( 27) 

B ( 22 ) 

A ( 83 ) 
6(5) 

B ( 66 ) 

A ( 46 ) 

B ( 32 ) 

B ( 36) 

B( 37 ) 

B ( 1 5 ) 

B( 177 ) 
A ( 1 53 ) 
A( 123) 

Ad) 

A ( 60 ) 

B( 51 ) 


JULIAN DATE OF TAKEOFF 

LARGEST OF THE RELATIVE ERRORS BETWEEN R-K AND LOW- 
ORDER INTEGRATION METHODS, EQ* (D4) 

LIST OF BCD BODY NAMES WHOSE POSITIONS ARE TO BE DE- 
TERMINED FROM TAPE DATA 

INITIAL ELEVATION ANGLE, USED WHEN IM0DE=4, SKETCH(B), 
DEGREES 

ELLIPSE DATA FOR PERTURBATING BODIES, READ FROM CAROS, 
12 PIECES OF DATA PER BODY 

ECCENTRICITY -1 

SEE TABLE II 

SQUARE ROOT OF ( ECCENTR ICITY SQUARED -1) 

SEE TABLE II 
SEE TABLE II 

NATURAL LOGARITHM OF EREF 
SEE TABLE II 

AEXIT ( NSTAGE ) / 100, NEWTONS/MB 
ECCENTRICITY CALCULATED WHEN IM0DE=3 
SEE TABLE II 

RATE OF PROPELLENT FLOW, KG/SEC 
FL0W1 I NSTAGE ) 

X, Y, Z COMPONENTS OF THRUST ACCELERATION, M/SEC**2 
DEFINED IN SUBROUTINE AERO, SET IN STDATA 
GEOPOTENTIAL, M 

GRAVITATIONAL CONSTANT, MU, OF THE SYSTEM, M**3/SEC**2 

SQUARE ROOT OF GK2M 

VALUE OF DELT FOR PREVIOUS STEP 

DEFINED IN SUBROUTINE ORDER 

SEE TABLE II 

INPUT IDENTIFICATION NUMBERS ASSOCIATED WITH EACH 
STAGE 

SEE TABLE II 

SET OF INDICES, SET IN STDATA 

NUMBER OF SETS OF ERROR DATA, SET IN ERRORZ FOR USE 
IN NBODY 
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INLOOK 


A ( 599 ) 

INPUT IDENTIFICATION NUMBER FOR INPUT AFTER FINDING C 
(LOOKX) = XLOOK 

KSUB 


B ( 19 ) 

INDEX OF RUNGE-KUTTA SUBINTERVALS 

fcAT 


A ( 33) 

INITIAL GEOCENTRIC LATITUDE, USED WHEN IM0DE=4, SKETCH 
(B), DEGREES 

LONG 


A I 34) 

INITIAL LONGITUDE RELATIVE TO GREENWICH, USED WHEN 
I MQ0E=4* SKETCH! B ) , DEGREES 

LOOKX 


A( 8 ) 

SEE TABLE II 

LOOKSW 


A ( 9 ) 

SEE TABLE II 

LSTAGE 


At 38) 

TOTAL NUMBER OF STAGES INTEGRATED BEFORE RETURNING TO 
THE MAIN PROGRAM 

MBODYS 


B (42) 

NUMBER OF PERTURBATING BODIES (NBODYS-1) 

HODOUT 


A (20) 

SEE TABLE II 

NBODYS 


B ( 41 ) 

TOTAL NUMBER OF BODIES, EXCLUDING THE VEHICLE 

NCASES 


A ( 600 ) 

SAVED VALUE OF NCASE 

NCASE 


cm 

CASE NUMBER, RAISED ONCE EACH TIME CONTROL PASSES 
THROUGH THE MAIN PROGRAM 

NEFMRS 

(8) 

B ( 1 85 ) 

DEFINED IN SUBROUTINE ORDER 

NEQ 


A ( 2 ) 

NUMBER OF EQUATIONS TO BE INTEGRATED, SET TO 8 IN 
STDATA 

NSAVE 


C(4) 

SEE TABLE II 

NSTAGE 


AC3) 

THE INDEX INDICATING THE PARTICULAR STAGE CURRENTLY 
BEING INTEGRATED 

NSTART 


B ( 24 ) 

INTERNAL CONTROL IN NBODY AND EQUATE 

OBLATJ 


AC 26 ) 

OBLATENESS COEFFICIENT OF SECOND HARMUNIC 

OBLATD 


A ( 27 ) 

OBLATENESS COEFFICIENT OF FOURTH HARMONIC 

OBLATH 


A ( 28 ) 

OBLATENESS COEFFICIENT OF THIRD HARMONIC 

OBLATN 


A C 40 ) 

SEE TABLE II 

OBLAT i 

13) 

B< 75) 

X, Y , Z COMPONENTS OF OBLATENESS ACCELERATION, M/SEC**2 

OLDDEL 


B ( 9 ) 

VALUE OF DELT FOR PREVIOUS GOOD STEP 

ORBELS 

16) 

BC 1 16 ) 

ARRAY OF OUTPUT VAR I ABLES , E I THER RECTANGULAR OR ORBIT 
ELEMENTS 

OUTPOT 


BC 399) 

CAUSES ABSENCE OF OUTPUT WHEN NONZERO 

P ( 31 


B( 84 ) 

DEFINED IN EQ* (B4> 


29 


PAR 43) 

B ( 60 ) 

DEFINED BY EQUATIONS IN SUBROUTINE THRUST 

PMAGN 


B ( 50 ) 

0EF1NED IN EQUATION FORM BY SUBROUTINE THRUST 

PNAME 

(30) 

A (287 ) 

PERMANENT LIST OF BODY NAMES MADE FROM PNAA LIST IN 
SUBROUTINE ORDER* ELIPS NAMES BEGIN AT PNAME! 21 ) 

PRESS 


B ( 33 ) 

ATMOSPHERIC PRESSURE* MB 

PSI 


BOO) 

PATH ANGLE* ANGLE BETWEEN PATH AND LOCAL HORIZONTAL* 
DEGREES 

PSIR 


B ( 398 ) 

RELATIVE PATH ANGLE, TAKEN RELATIVE TO A ROTATING 
ORIGIN BODY, DEG 

PUSH 


A ( 166) 

THRUST FORCE, NEWTONS 

PUSHO 


B09L) 

VACUUM THRUST FORCE, NEWTONS 

Q 


B( 59 ) 

DYNAMIC PRESSURE* NEWT0NS/M**2 

UMAX 


B ( 44) 

MAXIMUM VALUE OF Q DEVELOPED DURING A SINGLE TRAJEC- 
TORY ( SET TO ZERO WHEN CONTROL PASSES THROUGH SUB- 
ROUTINE EXTRA) 

QX (3) 


B( 78 ) 

X* Y« Z COMPONENTS OF PERTURBATIVE ACCELERATION DUE TO 
PERTURBATING BOOIES* M/SEC**2 

R (8) 


B( 102 ) 

DISTANCES OF ALL BODIES FROM OBJECT, IN ORDER OF 
bname LIST, m 

RADIAL 


BC81) 

RADIAL COMPONENT OF TOTAL PERTURBATIVE ACCELERATION, 
POSITIVE OUTWARD, M/SEC**2 

RAMC ( 5 ) 

B ( 393 ) 

RELATIVE ANGULAR MOMENTUM PER UNIT MASS COMPONENTS, 
TOTAL RELATIVE ANGULAR MOMENTUM PER UNIT MASS, AND 
ITS SQUARE, M**2/ SEC 

RATM 


A ( 22 ) 

RADIUS OF ATMOSPHERE, M 

RATMOS 


B( 23 ) 

SET EQUAL TO RATM WHEN ATMN EQUALS THE REFERENCE BODY 
NAME, BNAME ( I ) 

RATIO 


B ( 58 ) 

RATIO OF ADJACENT STEP SIZES, DELT 

R6 (3, 

8) 

B ( 193) 

X, Y , Z COMPONENTS OF DISTANCE FROM ALL BODIES TO THE 
OBJECT, M 

RBCRIT 

(8) 

B ( 145) 

LIST OF SPHERE-OF-INFLUENCE RADII OF ALL BODIES IN 
BNAME LIST, M 

RCRIT 

(30) 

A( 377 ) 

PERMANENT LIST OF SPHERE-OF-INFLUENCE RAOII CORRES- 
PONDING TO PNAME LIST OF BODY NAMES. RADII FROM ELIPS 
DATA BEGIN AT RCRITI21), M 

RE 


A ( 25 ) 

RADIUS OF EARTH EQUATOR, M 

RECALL 


C ( 5 ) 

SEE TABLE II 

REFER 

(30) 

A ( 317 ) 

LIST OF REFERENCE BODIES CORRESPONDING TO PNAME LIST, 
REFERENCE BODIES FROM ELIPS OATA BEGIN AT REFER! 21 ) 

RESQRD 


B( 7 ) 

SQUARE OF RE 

RETURN 


B ( 400 ) 

CAUSES CONTROL NOT TO RETURN TO MAIN PROG. IF NONZERO 
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mu hi in i 



I 


REVS 

REVOLV 

RMASS1 (10) 

ROTATE 

RSQRO 

SIGNAL 

SIMP1 (10) 

SIMP 

SINALF 

SINBET 

SINTRU 

SINCE 

SINV 

SPACES 

SPD 

SQRDKi 

SQRDK 

STEPNX 

STEPS 

STEPGO 

STEPNO 

SWLOOK 

TABLT 

TABLE (200) 

TAPE3 
TB (10) 


A ( 48 ) 
8 ( 21 ) 

A ( 73 ) 

A ( 39 ) 

B ( 45 ) 

B ( 31 ) 

A ( 93 ) 

B ( 2 ) 
8(46) 

B ( 47 ) 

B ( 52 ) 

B ( 54 ) 
6(56) 

B( 16) 

A ( 44) 

A ( 47 ) 

B (35) 

A( 16 ) 

A ( 1 7 ) 

A t 4 1 ) 

A ( 42 ) 

A( 10) 

B ( 20 ) 

C ( 191 1 ) 

C ( 2 ) 

A ( 63 ) 


REVOLUTION COUNTER, USED ONLY FOR OUTPUT 

ROTATION RATE OF REFERENCE BODY WHEN ATMN=BNAME( 1) , 
RAD/SEC 

INITIAL MASSES FOR AT MOST 10 STAGES, KG 
ROTATION RATE OF A REFERENCE BODY, RAD/SEC 
RADIUS SQUARED OF OBJECT TO ORIGIN, M**2 
SEE TABLE II 

SPECIFIC IMPULSES FOR AT MOST 10 STAGES, SEC 
S I MP 1 ( NST AGE ) 

SINE OF ALPHA 

SINE OF BETA 

SINE OF TRU 

SINE OF INCLINATION 

SINE OF THE ARGUMENT OF LATITUDE 

NUMBER OF EQUAL TIME UNITS UNTIL NEXT OUTPUT 

SECONDS PER DAY, SET IN STDATA, SEC/DAY 

GRAVITATIONAL CONSTANT OF THE SUN, AU**3/DAY**2 

GRAVITATIONAL CONSTANT OF THE SUN, M**3/SEC**2 

SEE TABLE II 

SEE TABLE II 

COUNT OF SUCCESSFUL INTEGRATION STEPS 

COUNT OF UNSUCCESSFUL INTEGRATION STEPS (THOSE WHICH 
DO NOT PASS ERROR CONTROL TEST) 

SEE TABLE II 

TIME MEASURED RELATIVE TO THE JULIAN DATE OF TAKEOFF, 
DAYS 

ARRAY OF INPUT PARAMETERS AND THEIR COMMON STORE LO- 
CATIONS 

SEE TABLE II 

FLIGHT TIMES FOR AT MOST 10 STAGES, SEC 
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TDATA (6,3,7) 

B ( 265 ) 

COEFFICIENTS FROM EPHEMERIOES TAPE TO BE USED IN DE- 
TERMINING POSITIONS AND POSSIBLY VELOCITIES OF PER- 
T UR BATING BODIES, ONE SET FOR EACH OF 7 BODIES 

TDEL (7) 

B( 170) 

ONE-HALF OF TIME SPACING BETWEEN TWO ADJACENT ENTRIES 
OF LIKE BODY NAME ON EPHEMERIOES TAPE, READ FROM 
TAPE FOR EACH BODY 

TFILE 

A ( 6 ) 

SEE TABLE II 

TIM (7) 

B ( 163 ) 

TIME FOR SET OF EPHEMERIS DATA, READ FROM EPHEMER IDES 
TAPE, ONE FOR EACH BODY 

TKICK 

A( 15) 

INITIAL STEP SIZE OF A TRAJECTORY TO BE COMPUTED IN 
CLOSED-FORM, FOR USE WHEN IM0DE=4, WHICH FACILITATES 
STARTING OF SOME TYPES OF TRAJECTORIES 

TM 

BI34) 

ATMOSPHERIC TEMPERATURE TIMES THE RATIO OF MOLECULAR 
TO ACTUAL MOLECULAR WEIGHT, DEGREES KELVIN 

TMAX 

6(4) 

SEE TABLE II 

THIN 

A( 18) 

SEE TABLE II 

TOFFT 

A ( 24 ) 

FRACTIONAL PART OF JULIAN OATE OF TAKEOFF, DAYS 

TRSFFR 

B C 8 ) 

SEE TABLE II 

TRU 

B ( 40 ) 

TRUE ANOMALY, RAD 

TTEST 

A ( 54 ) 

SEE TABLE II 

ITOL 

A ( 45 ) 

TIME TOLERANCE WITHIN WHICH PROBLEM TIME MINUS TMAX 
MUST LIE TO END STAGE 

U 

A ( 59 ) 

ECCENTRIC ANOMALY, RAD 

V 

B ( 95 ) 

VELOCITY OF OBJECT RELATIVE TO THE ORIGIN, M/SEC 

VATM (3) 

B ( 97 ) 

X, Y, Z COMPONENTS OF THE RELATIVE VELOCITY, VQ, M/SEC 

VEFM (3,8) 

B ( 241 ) 

X, Y, Z COMPONENTS OF OBJECT VELOCITY RELATIVE TO ALL 
BODIES, M/SEC 

VEL 

A ( 37 ) 

INITIAL RELATIVE VELOCITY, USED WHEN IMODE=A, SKETCH 
(B), M/SEC 

VMACH 

B ( 38 ) 

MACH NUMBER OF OBJECT 

VQ 

B ( 100 ) 

VELOCITY OF OBJECT RELATIVE TO ATMOSPHERE, M/SEC 

VQSQRD 

8(101) 

SQUARE OF VQ, M**2/SEC**2 

VSQRD 

B ( 96 ) 

SQUARE OF V, M**2/SEC**2 

VX 

B ( 92 ) 

X COMPONENT OF VELOCITY, M/SEC 

VY 

B ( 93 ) 

Y COMPONENT OF VELOCITY, M/SEC 

VZ 

B ( 94 ) 

Z COMPONENT OF VELOCITY, M/SEC 

X (100) 

B ( 401 ) 

WORKING SET OF INTEGRATION VARIABLES 

XDOT (100) 

B ( 501 ) 

TIME DERIVATIVES OF THE SET X 
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X I FT (3) 

B ( 72 ) 

X » Y * Z COMPONENTS OF LIFT ACCELERATION, M/SEC**2 

XINC (100) 

B ( 60 1 ) 

INCREMENTS OF THE INTEGRATION VARIABLES PER STEP 

XL00K 


A ( 12) 

SEE TABLE II 

XP (3, 

8) 

B ( 2 17 ) 

X, Y, Z COMPONENTS OF PERTURBATING BODY POSITIONS RELA- 
TIVE TO ORIGIN 

XPRIM 

(100,2) 

C ( 7 1 1 ) 

TWO 100-ELEMENT SETS, THE FIRST SET CONTAINS VALUES 
OF THE INTEGRATION VARIABLES AT THE PREVIOUS GOOD 
STEP, THE SECOND SET IS UNDER THE INTEGRATION PROCESS, 
SEE TABLE V 

XPRIMB 

(100,2) 

C ( 9 1 1 ) 

LEAST SIGNIFICANT HALF OF DOUBLE PRECISION INTEGRA- 
TION VARIABLES XPRIM 

XTOL 


A( 11) 

TOLERANCE ON THE DI SC IM INAT I ON C< LOOKX)-XLOOK TO BE 
SATISFIED 

XWHOLE 

(6) 

B ( 1 10 ) 

RECTANGULAR COORDINATES AND VELOCITIES, SET ASIDE FOR 
USE IN ORIGIN TRANSLATIONS 

ZN 


BC43) 

MEAN ANGULAR MOTION OF OBJECT, RAD/SEC 

ZORMAL 


B(83) 

Z COMPONENT OF TOTAL PERTURBATIVE ACCELERATION, 
M/SEC**2 
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APPEMDIX F 


LEWIS RESEARCH CENTER EPHEMERIS 
General Description 

The ephemeris data initially available on magnetic tape were from the Themis 
code prepared by the Livermore Laboratory, evidently from U. S. Naval Observatory 
data* Later, an ephemeris was obtained from the Jet Propulsion Laboratory assem- 
bled as a joint project of the Jet Propulsion Laboratory and the Space Technology 
Laboratory. These data are given relative to the mean vernal equinox and equator 
of 1950.0 and are tabulated with ephemeris time as the argument. 

An ephemeris was desired for certain uses in connection with the IBM 7090 
computer that would be shorter than the original ephemeris tapes mentioned and 
would be as accurate as possible consistent with the length. A short investiga- 
tion of the various possibilities led to adoption of fitted equations. In par- 
ticular, fifth-order polynomials were simultaneously fitted to the position and 
velocities of a body at three points. This procedure provides continuity of po- 
sition and velocity from one fit to the next, because the exterior points are 
common to adjacent fits. Polynomials were selected rather than another type of 
function, because they are easy to evaluate. Three separate polynomials are used 
for the x, y, and z coordinates, respectively. 


Procedure Used to Fit Data 


The process of computing the fitting equations is as follows; 

(l) A group of 50 sets of the components of planetary position was read 
into the machine memory for a single planet together with differences as they 
existed on the original magnetic tape. The differences were verified by compu- 
tation (in double precision because some data required it); and any errors were 
investigated, corrected, and verified. Published ephemeris data were adequate 
to correct all errors found. 


(2) The components of velocity v x , Vy, and v z were computed and stored 
in the memory for each of the 50 positions by means of a numerical differentia- 
tion formula using ninth differences; namely. 


X = (Ti - T_ x ) 


Al_i + AI + ]_ 

2 


ACII.! + AIII +1 AV_ ± + AV +1 
12 + 60 


AVTI.-l + AVII +1 AIX.-l + Aix +1 
" 280 + 1260 “ ^ F1 

(See ref. 11, PP* 42 and 99 for notation. ) Double-precision arithmetic was used 
for differences, but velocities were tabulated with single precision. 
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(3) Coefficients C , D, E, and F in the fifth-order polynomial 

X = Xq + S^(T - T 0 ) + C(T - T q ) 2 + D(T - T Q ) 3 + E(T - T Q ) 4 + F(T - T Q ) 5 (F2) 

and its derivative 

X = + 2C(T - T 0 ) + 3D(T - T Q ) 2 + 4E(T - T Q ) 3 + 5F(T - T 0 ) 4 (F3) 

were found to fit a first point (which was far enough from the beginning point 
to have all differences computed) and two equally spaced points for each compo- 
nent of position and velocity* (The initial spacing is not important, as will be 
seen later. ) Spacing is defined as the number of original data points fitted by 
one equation. Single-precision arithmetic was used. 

(4) The coefficients C, D, E, and E in step (3) were then used in equa- 
tions (F2) and (E3) to calculate components of all positions and velocities given 
in the original data and lying within the interval fitted. These values were 
checked with the original data. Radius R and velocity V were computed at the 
times tabulated in the original data. If any component of the position differed 
from the original data by more than RX1CT 7 or if any velocity differed from the 
original by more than 7X10-6, the fit was considered unsatisfactory. 

(5) If the fit was considered unsatisfactory, this fact was recorded, and 
the spacing was reduced by two data points. Steps 2 to 4 were then repeated. If 
the fit was considered satisfactory, this fact was recorded, and the spacing was 
increased by two spaces. Steps 2 to 4 were repeated. The largest satisfactory 
fit was identified when a certain spacing was satisfactory and the next larger 
fit was not satisfactory. 

(6) The coefficients that corresponded to the largest satisfactory fit were 
recorded on tape in binary mode as follows: 


Word 

Data 

Mode 

Definitions and/or units 

number 




1 

Planet name 

BCD 

Six characters (first six) 

2 

Julian date 

Floating point 

Date of midpoint of fit, Julian date 

3 

Delta T 



Number of days on each side of midpoint 

4 

F 



a AU/ day5 

5 

Ex 



a AU/day4 

6 

D x 



a-AU/day3 

7 

c x 



a AU/day2 

8 

* 



a AU/day 

9 

X 



a AU 

10 

Fy 



a AU/day 5 

11 

4 



a AU/day 4 

12 

% 



a AU/day3 

13 



a AU/day 2 

14 

y 



a AU/ day 

15 

y 



a AU 

16 

17 

l\ 



a AU/ day 5 
a AU/ day 4 

18 

D z 



a AU/ day 5 

19 

0 Z 



a AU/day2 

20 

i 



a AU/ day 

21 

z 



a AU 


a Except for Moon data, which are in Earth radii and days. 
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(7) As soon as a set of coefficients -was selected for an interval, addi- 
tional data were read from the source ephemeris tape and used to replace the 
points already fitted (except the last point). These data were processed as de- 
scribed in steps 1 and 2 so that the next 50 points were ready to be fitted. 
Steps 3 to 6 were then used to find the next set of coefficients, and steps 1 
to 6 were repeated until all data for all planets were fitted. 


Data Treated 

The preceding process was applied to all data available at the time. For 
the Moon, the technique usually led to the use of every point in the fitted in- 
terval (i. e. , only three points were fitted). Thus, a check of accuracy was not 
available. The error in the attempt to fit the next greater interval (five 
points) was not excessive, however, and it is .judged that the accuracy obtained 
from these fits is about equal to that held on the other bodies. 


Merged Ephemeris Tape 

Once all the positions and velocities of all the bodies then available were 
fitted, the coefficients were merged in order of the starting date of each fit. 
The resulting tape was written in binary mode with 12 sets of fits per record. 

The detail of this record is as follows: 


Set 1 


Set 2 


1st word: 
2nd word: 
3rd word: 
4th word: 


23rd word: 


PORTRAIT compatible 

file number, fixed point in decrement 

planet name, code in BCD, first six characters 

Julian date, floating point 

etc. , according to list in paragraph 6 

21 words 

z 


24th word: 
25th word: 


planet name, code in BCD, first six characters 
Julian date, floating point 


.44th word: 


z 


Successive sets follow one another with a total of 12 sets. 

Set 12 f 234th word: planet name 

(last set) •S 235th word: Julian date, floating point 

(.254th word: z 

End-of-record gap 


One record contains 254 words, the first is for FORTRAN compatibility, the second 
is a file number used for identification in the system. It is a fixed point 2. 
The third is the beginning of the first set of data, and 12 sets follow each with 
21 words. The last word is the 254th word (counting the FORTRAN compatible word) 
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followed by an end-of-record gap* The remaining records are compiled in the same 
manner with an end-of-file recorded as a terminating mark* 

Because of the merging operation, all bodies are given in one list in a 
random order according to the starting date of the interval. The starting date 
is the Julian day (word 2) minus the half interval (word 3) (see procedure, par- 
agraph 6). The entire ephemeris occupies about one-seventh reel of tape. A 
summary of data is given in table VII* 
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APPENDIX G 


INPUT-DATA REQUIREMENTS 

The procedure needed to run actual problems with the aid of this routine is 
described herein. It is intended to permit the user with a specific problem in 
mind to make a complete list of data required and to select desirable operating 
alternatives from those available. The details of this procedure are contained 
in the following instructions! 

(1) Provision has been made for two types of ephemeris data to specify the 
locations of celestial bodies that perturb the vehicle. They are ellipse data 
and ephemeris-tape data. If the problem does not involve perturbing bodies (ex- 
cept a reference body) or if elliptic data are used for all the perturbing 
bodies, skip to instruction 5. 

(2) If the perturbing-body data are to be taken from an ephemeris tape, list 
the names of the ephemerides and Julian dates to be covered along with the fol- 
lowing auxiliary information: 

1st card: $DATA = 300, $TABLE, 2 = TAPE 3, 17 = ELIST, 29 = TBEGIN, 

30 = TEND/ 

Other cards : TAPE 3=0 

TBEGIN = ephemeris beginning Julian date 
TEND = ephemeris ending Julian date 

ELIST = (names of perturbing bodies in "ALF" format, see 
example in text) 

The ephemerides of all planets except Earth bear the name of the planet. The 
ephemeris giving the distance from Earth to the Sun is called "sun," as is 
astronomical practice. 

(3) If successive files on the ephemeris tape are to be made, punch the 
corresponding sets as follows: 

$DATA = 300, TAPE 3=0, TBEGIN = , TEND = , ELIST = 

As many similar sets as are needed may be appended. 

(4) If ellipse data are to be loaded from cards, they are prepared later 
under instruction 11. 

(5) On the first execution after loading the routine, the common area is 
cleared whether an ephemeris tape is constructed or not. It is now necessary to 
load a table of variable names. Once loaded, this table will not be cleared 
again (except if the control variable TAPE 3 is set equal to zero). These names 
are for use on the input cards. If a different name is desirable for any 
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variable, it may be- changed in the table and where it appears on the input card 
(ref. 7). The cards are: 

$ DAT A=l, Ip TABLE, 33=DTOFFJ, 34=TOFFT, 7 UPTIME, 716=X, 717=Y, 718=Z, 
713=VX, 714=VY, 715=VZ, ll.=IMODE, 713=E, 714=0MEGA, 715=N0DES, 
716=INCL, 717=MA, 718=P, 43=LAT, 44=L0NG, 45=AZI, 46=ELEY, 14=ALT, 
47=YEL, 16=TFILE, 28=TMIN, 153=B0DYCD, 177=ELIPS, 30.=M0D0IJT, 
27=STEPS, 29=DELMAX, 26=STEPMX, 23=EREF, 24=ERLIMT, 4.=NSAVE, 
5=RECALL, 3=CLEAR, 18.=L00KX, 22=XL00K, 19.=LOOKS¥, 20=SWL00K, 

609 . =INL00K, 15=END, 31=ATMN, 32=RATM, 49=R0TATE, 417=C0EFN, 
163.=ICC, 60=BETA, 50=0BLATN, 73=TB, 93=FL0W, 103=SIMP, 123=AREA, 
143=DELT, 83=RMASS, 113=AEXIT, 133.=IDENT, 48.=LSTAGE, 25=TKICK / 

(6) The initial position and velocity of the vehicle may be given in any one 
of the three coordinate systems. If the initial data are given in orbit ele- 
ments, skip to instruction (8). If the initial data are given in rectangular co- 
ordinates, skip to instruction (7). If the initial data are given in Earth- 
centered spherical coordinates, the following variables should be punched; 

LAT « latitude, deg, positive north of equator 

LONG = longitude, relative to Greenwich, deg 

ALT = altitude above sea level, m 

AZI = azimuth angle, east from north, deg 

ELEV* - elevation angle, horizontal to path, deg 

VEL = initial relative velocity, m/sec 

TKICK = size of initial vertical, nondrag step to facilitate starting, 
sec 

If the Earth is assumed to be rotating but aerodynamic forces are not 
to be considered, set 

ROTATE - Earth rotation rate, 7.29211585X10“^ radian/sec 

If integration in rectangular coordinates is desired set 
IMODE = 4 

or else if integration in orbit elements is desired set 
IMODE = -4 

Skip to instruction (9). 

(7) If the initial data are in rectangular coordinates, set the following 
variables : 

X = x-component of position in x,y, z coordinate system, m 
Y y-component of position in x,y, z coordinate system, m 
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Z = z-component of position in x,y,z coordinate system, m 

VX = x-component of velocity in x,y, z coordinate system, m/sec 

VY = y-component of velocity in x,y,z coordinate system, m/sec 

VZ = z-component of velocity in x,y,z coordinate system, m/sec 

If integration in rectangular coordinates is desired set 
IMODE = 2 

or else, if integration in orbit elements is desired set 
IMODE = -2 

Skip to instruction (9). 

(8) If the initial data are in orbit elements, set the following variables: 

E = eccentricity 

OMEGA = argument of pericenter, radians 

NODES = longitude of ascending node (to mean vernal equinox of 1950.0), 
radians 

INCL s= orbit inclination to mean equator of 1950.0, radians 
MA = mean anomaly, radians 
P = semilatus rectum, m 

If integration in orbit elements is desired set 
IMODE = 1 

or else, if integration in rectangular coordinates is desired set 
IMODE = -1 

(9) To specify takeoff time, set the following variables: 

DTOFFJ = Julian day number 
TOFFT = fraction of day 

TIME = time from previously set Julian date, sec 

Takeoff occurs at the instant (ephemeris time) corresponding to the sum of the 
last three quantities. If a specific date or time is not required, these vari- 
ables may be skipped. In that case, the SUBROUTINE STDATA sets DTOEFJ to 
2440 000. 

(10) To specify the origin and any perturbing bodies, list them as BODYCD = 
(list of body names in U ALF" format, see text example). The first body in the 
list is taken to be the reference body. The distances between the bodies in 
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this list must he computable from either ellipse data (instruction (ll) ) or 
ephemeris-tape data (instruction (2)). There may be no more than eight names in 
the list* Also, if the ephemeris tape is being used, the correct file must be 
found on it. For this purpose, set TFILE = desired ephemeris-tape file* The 
ephemeris files were numbered in sequence when written in instruction (2). If 
TFILE is not given, it will be set equal to 1.0 by the SUBROUTINE STDATA. 

(11) For each body whose path is represented by an ellipse, a 12-element 
set of data must be loaded. A 12-element set consists of: 

1. Body name in "ALF M format (maximum of six characters) 

2* Reference body name in "ALF" format (maximum of six characters) 

3* Mass of body, sun mass units 

4* Radius of sphere of influence, m 

5. Semilatus rectum, AU 

6. Eccentricity 

7* Argument of pericenter, radians 

8. Longitude of ascending node (to mean vernal equinox of 1950.0), 

radians 

9. Orbit inclination (to mean equator of 1950.0), radians 

10. Julian day at perihelion 

11. Fraction of day at perihelion 

12. Period, mean solar days 

It is convenient to punch a 12-element set in sequence and to separate the ele- 
ments by commas on as many cards as are required. Several sets may then be 
loaded consecutively. The order of the sets is immaterial. Ellipse data, if 
present, take precedence over ephemeris-tape data. The sets are loaded consecu- 
tively, in any order, as follows: 

SLIPS = set 1, set 2, set 3, . . set nj n < 10 (see example in appen- 
dix I) 

(12) If oblateness effects of the Earth are to be included, set 
0BLATN = (ALF5) EARTH 

(13) Provision has been made to fly multistage vehicles with up to 

10 stages. At least one stage must be loaded. There are eight parameters for 
each stage with provision for input- controlled modifications of other variables. 
The 10 values of each parameter are stored in an array corresponding to the 
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10 stages. Input cards are as follows : 

GIB = burning time for 1st stage, 2nd stage, etc. , sec 

FLOW = propellant flow rate for 1st stage, 2nd stage, etc. , kg/sec 

SIMP = vacuum specific impulse of 1st stage, 2nd stage, etc. , sec 

AREA — aerodynamic reference area of 1st stage, 2nd stage, etc. , m^ 

AEXTT » engine exit area for 1st stage, 2nd stage, etc. , vc? 

EMASS = initial mass or jettison mass for 1st stage, 2nd stage, etc. , 
kg 

DELT = initial integration step size for 1st stage, 2nd stage, etc., 
sec 

IDEM 1 = input identification number 1st stage, 2nd stage, etc. 

GIB must be loaded for as many stages as are to be flown. Others may be omitted 
if zero is appropriate. If EMASS(i) is not positive, the i th stage begins with 
the final mass of the previous stage reduced by the fixed amount EMASS(i). In 
the case of PELT, zero will result in use of GIB/lOO. IDEKT of a nonzero value 
will cause any data cards of that identification number to be read in after the 
stage is set up and before integration begins. Gdris permits the user to make 
almost any change desired. GIhe order of data cards is discussed in instruc- 
tion (24). 

(14) G?he thrust orientation must be specified by setting 
BETA » angle p, deg (see sketch (a) (p. 4)) 

COEFR (I) = angle-of-attack schedule, a = a(t) (see instruction (16)) 

ICC = fixed-point integer (see instruction (16)) 

For the special case of tangential thrust, none of the last three variables need 
be set. 

(15) If aerodynamic forces are present, set in addition to AREA in instruc- 
tion (13) : 

AGMKT = name of body that has atmosphere, in "ALF" format, (Earth) 

RAGM = radius above which atmospheric forces are not to he considered, 
m 

ROTATE = atmospheric-rotation rate, radians/sec ( 7. 29211585X10"^ for 
Earth) 

BETA = angle (3, deg (see sketch (a)) 


J 
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COEFE (i) = angle-of-attack schedule, a, = a(t), Cx/sin a, Cj^Q/ an( ^ 
C^ curves (see instruction (16)) 

ICC = fixed-point integers (see instruction (16)) 


(16) If neither thrust nor aerodynamic forces are present, skip to instruc- 
tion (18). The relations a(t), C T /sin cl, C^ q, and C^ . / C^ are assumed to he 


D,i 7 L 


quadratic functions that involve coefficients, which are located in the COEFE(J) 
array. The arrangement of these coefficients is best explained by an example. 
Suppose the function a(t) is as follows: 


a ll 

a 21 


'‘izt 

+ a 13 t 2 

(t-L < t < t 2 ) 

'•2Z t 

+ a 23 t 2 

( *fc 2 ^ *fc ^ 1j 3 ) 

l 5Z t 

+ a 33 t2 

(t 3 < t < t 4 ) 


etc. 


etc. 


The coefficients a.? ± should then be loaded into the COEFE(J) array as: 
COEPN(J) = ti , &2 -| > a 22’ a 23’^3> a ‘'51 )S 'Z2’ a '3Z’^4:’ 


• ♦ • 




Furthermore, additional sets of coefficients for the other functions may simply 
be added to the COEFN(J) array, which results in a string of sets of coeffi- 
cients^ and can be represented, for example, as: 

COEFE(J) = a coefficients, C^/sin a coefficients, Cj^q coefficients, etc. 

= t 1 ,a 11 ^a 12 ^ . • • etc * 


The starting point in the COEFN(J) array of each function must also be loaded to 
identify the correct region of coefficients* To this end, the following array 
must also be loaded: 


ICC(l) = fixed-point value of 


J where a coefficients begin 


ICC (2) ss fixed-point value of J where 
ICC ( 3 ) = fixed-point value of J where 


C^/sin a coefficients begin 
c D,i/^L coefficients begin 


ICC(4) =s fixed-point value of J where C-q^q coefficients begin 

For this purpose, all values in the COEFN(J) array are called coefficients (i.e. , 
the t 1 s and the s are coefficients). The sequence of the sets is arbi- 
trary, since changing the sequence requires only a change in the ICC(l) array. 
(See appendix I for Example II, the lunar orbiting probe.) 
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•(17) The size of the integration steps is determined primarily hy the error 
control variables* These are loaded asi 

EEEF 5 = error reference value j 15 in appendix D 

EELIMT = maximum value of 5 that is acceptable on any particular step 

EEEF is always treated as a positive number^ however, if it is loaded with a 
minus sign., this will cause error information to be printed at the completion of 
the problem* If no error control data is loaded, SUHRO'uTlSE STDATA will set 
EEEF = 2X10-6, EKLIMT = 3X10“ 6 . 

(18) The output control offers a choice on the frequency of output data as 
follows ; 

If M0D0UT = 1, output will occur every n^ step (n = STEPS) until 
t s* TMTTT, and then M0D0UT is set equal to 2 by the program 

If M0D0UT = 2, output occurs at equal time intervals of DEIMAX until 
t = TMAX 

If M0D0UT = 3, output occurs at equal time intervals of DEIMAX until 
t = BUN, then M0D0UT is set equal to 4 by the program 

If M0D0UT = 4, output occurs every n^* 1 step (n = STEPS) until 
t = TMAX 

STEPMX = maximum step limit before problem is completed 
DEIM AX = time interval between outputs 
STEPS = number of steps between outputs 
3MHT = t-ime when M0D0UT changes 

Note that output control may, at times, strongly influence the integration step 
size especially if M0DOUT is 2 or 3 and DEIMAX is small. STDATA will put 
M0D0UT = 4 and STEPS = 1. 

Note that TMAX «* time at start of a stage, plus the stage time, TB(NSTAGE), and 
is computed internally. 

(19) Provision has been made to interrupt the integration procedure when an 
arbitrary value of an arbitrary parameter is attained. By interrupt it is meant 
that an output will occur at this point, input is permissible, and a decision is 
made whether to continue the stage, terminate the stage, or terminate the flight. 
Skip to instruction (20) if this facility is not desired. To cause an interrupt, 
set 

L00KX = COMMON C location of arbitrary parameter 
XL00K = value of C(LOOKX) where an interrupt is desired 
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INLOOK = input identification number for interrupt 


END = a negative number if flight should be terminated, zero if stage 
should continue, or a positive number if stage should be terminated 

If the interrupt is not desired the first time C(LOOKX) ~ XLOOK, set 

LOOKSW as COMMON C location of a second arbitrary parameter 

SWLOOK « value of C (LOOKSW), which must be equaled or exceeded before 
an interrupt may occur (interrupt occurs if C(LOOKX) 2 = XLOOK and 
C (LOOKSW) > SWLOOK) 

Typically, time may be the second arbitrary parameter } thus, STRATA sets 
LOOKSW - 711, the COMMON location of time. INLOOK of a nonzero value will cause 
any data cards of that identification number to be read-in prior to the interro- 
gation of END. The order of the cards is discussed in instruction (24). 

(20) Provision has been made to save a block of initial conditions and pro- 
gram control parameters prior to the integration of the n^* 1 stage. This allows 
the flight to be flown again from the n th stage onward with prescribed altera- 
tions. Skip to instruction (2l) if this facility is not desired. To save the 
program control variable array. A, and the integration variable array, 

XPRIM + XPHIMB, just prior to integration of the n^* 1 stage, set 

NSAVE » the number of the n^* 1 stage 

The saved data, stored in the D array, will be returned to the A and 
XPRIM + XPRIMB arrays after the flight is completed if 

RECALL = any nonzero number 

It is intended that changes in the succeeding flight will be made at the main 
input station (^DATA=l). NSAVE and RECALL are not contained in the array A 
and are therefore unaffected by the save-recall sequence. The correct sequence 
of these controls is not always simple and an understanding of the main program 
and input stationing is quite desirable. 

(21) If the standard set of data contained in the StBROUTINE STRATA is not 
desired, set 

CLEAR any nonzero number 

It is intended that this control shall be set nonzero by the $DATA » 99 input 
station at the beginning of the main program. It is not affected by the save- 
recall sequencing explained in instruction (20). 

(22) If the number of stages to be flown is not equal to the number of con- 
secutive nonzero flight times, TB, set 

LSTAGE « number of last stage to be flown 
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(23) When a transfer of origin occurs, provision has been made to read input 
into the program* This is done with the aid of $DATA ** 101, followed by the 
data statements desired* 

(2 4) The sequencing of the input cards is not always simple and no rigid 
rules may be written down* Inspection of the program may be necessary to answer 
some questions* however, in general, the first input cards belong to the 
$DATA ss 300 group if an ephemerides tape is required* This group is followed 
by the $DATA « 1 group, which consists of the main input for a single flight. 
Following this are the in-flight input cards, if any, which may be any combina- 
tion of $r>ATA 5 = 101, $DATA - IHLOO K, or $DATA » XDEHT (TTSTAGE) groups* The 
order of these groups of cards matches the order of the time sequence of events 
in the flight itself* For multiple flights, sets of the above groups may be 
added in tandem* It is usually desirable in this case, however, to read all the 
$DATA =* 300 sets at the same time (as in instruction (3)) to avoid excessive 
tape handling* 

(25) Following is an input check list that may be helpful at execution time; 


INPUT CHECK LIST a 


1 Takeoff time 


Position and velocity 

1 Reference and perturbing bodies 



(completely fill in one and only one block) 

— 




— 




— 

BODYCD a 



Rectangular 

Orbit elements 

Spherical 

— 









Tape 

Elliptic 

DTOFFJ 

= 

X a 


E 

— 

LAT = 

b TAPE 3=0 

ELIPS a 

TOFFT = 


Y = 


OMEGA = 

LONG = 

b TBEGIN = 


TIME = 


Z « 


NODES = 

AZI a 

b TEND = 




VX a 


IN CL = 

EL EV a 

b ELIST = 




VY a 


MA = 

ALT a 

TFILE a 




VZ = 


P 


VEL a 





IMODE 

= 2 

IMODE = 1 

TKICK 









IMODE = 4 



Output 

Error 

Restart 

Parameter 

Atmosphere and 

Oblateness - 

Stage 

control 

control 

feature 

search 

coefficients 

rotation 

data 

TMIN = 

EREF a 

NSAVE - 

LOOKX = 

ATMN = 

OBLATN = 

TB = 

MODOUT ® 

ERLIMT = 

RECALL - 

XLOOK a 

RATM = 

ROTATE = 

FLOW a 

STEPS = 



CLEAR = 

LOOKSW = 

COEFN a 


SIMP = 

DELMAX » 





SWLOOK a 

ICC = 


AREA = 

STEPMX = 





INLOOK = 

BETA = 


DELT = 






END = 



RMASS = 









AEXIT = 









IDENT = 









LSTACE a 


a The following standard data are loaded by SUBROUTINE STDATA: 


DT0FFJ « 2440 000.0 
IMODE = 1 

BODY CD ( 1 ) a (ALF5) EARTH 
RMASS(l) « 1.0 

hAt input 300, setting TAPE 3 


MODOUT « 4 
STEPS « 1.0 
STEPMX » 100.0 
LOOKSW = 711 


EREF a 1x10“ 6 
ERLIMT « 3X10" 6 
TFILE « 1.0 


0 is necessary to make an ephemeris tape. 


i : 
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APPENDIX H 


PROGRAM LISTING 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


THIS MAIN PROGRAM IS THE SUPERSTRUCTURE ABOVE ALL SUBPROGRAMS. 

SUBROUTINE TAPE CLEARS COMMON 1 THRU 4000 AND MAY CONSTRUCT AN 
EPHEMER IS TAPE, ALSO, IT ALWAYS SETS T APES =0. SUBROUTINE STOAT A 
LOADS A STANDARD SET OF DATA. IF RECALL DOES NOT EQUAL ZERO , A 
PREVIOUSLY SAVED SET OF DATA ( FROM STAGE I IS MOVED TO THE INITIAL 
DATA LOCATION. THE MAIN INPUT STATION IS STATEMENT 8IINPUT1) 

WHERE THE VEHICLE DATA FOR ALL STAGES MAY BE LOADED. SUBROUTINE ORDER IS 
CALLED TO ORDER THE LIST OF BODIES, DETERMINE THE GRAVITATIONAL CONSTANT, 
ORIGIN ROTATION RATE, ATMOSPHERIC RADIUS, RELOCATE ELLIPTIC EPHEMER IS DATA 
AND POSITION THE EPHEMER I S TAPE. 

COMMON C 


DIMENSION A1600I, B(700), CC4000), 
1 TB(IO), DC1IOO) 


EQUIVALENCE 
1 1 A , C 

2ID ,C 

3INCASES, A 
MTB , A 


( ID), IB ,C 
(2111) ) > ( LSTAGE , A 
( 600) ) , (NiTAGE.A 
( 63) ) , (TAPE3 ,C 


( 1111D , (CLEAR ,C 
( 38) I , (NCASE ,C 
( 3) ) , (RECALL, C 
( 2) ) , (TABLE ,C 


( 3D, 

( D), 

I 5) ) , 
( 1911) ) 


1 CALL INPUT (99, C, TABLE) 

IF (TAPE3) 3,2,3 

2 CALL TAPE 

3 NCASE = NCASE + 1 

WRITE OUTPUT TAPE 6, 12, NCASE 

12 FORMAT! 12H1CASE NUMBER 1 3, 1H. ) 

IF (CLEAR) 5,4,5 

A CALL STDATA 

5 IF (RECALL) 6,8,6 

6 DO 7 J=l,1100 

7 A ( J ) = D(J) 

WRITE OUTPUT TAPE 6, 16 , NSTAGE , NCASES 

16 FORMAT ( 33H RECALLED INITIAL DATA FROM STAGE 12 , 8H OF CASEI4,1H.) 

8 CALL INPUT ( 1,C, TABLE) 

IF (SENSE Switch 6) 13,14 

13 WRITE OUTPUT TAPE 6,15 

15 FORMAT! 19HQEXIT VIA SENSE SW6> 

CALL EXIT 

14 IF (LSTAGE) 11,9,11 

9 DO 10 LSTAGh=l,10 

IF ( TB { L ST AGE + 1 ) ) 10,11,10 

10 CONTINUE 
LSTAGE = 10 

11 CALL ORDER 

17 CALL STAGE 
GO TO 1 
END 


SUBROUTINE AERO 

C SUBROUTINE AERO COMPUTES THE LIFT AND DRAG ACCELERATIONS. AS IN SUBROUT- 

C INE THRUST, THESE VECTORS ARE REFERENCED TO THE RELATIVE WIND VELOCITY. 

C COEFFICIENTS OF LIFT, INDUCED DRAG, AND DRAG AT ZERO ANGLE OF ATTACK ARE 

C ASSUMED TO BE FUNCTIONS OF MACH NUMBER AND ANGLE OF ATTACK. TABLES OF 

C CD I/CL**2» CL/SIN( ALPHA) , AND CDO ARE ASSUMED AS FITTED QUADRATIC EQUAT- 

C IONS IN THE COEFN ARRAY. GASFAC IS THE SQRTF ( SPEC I F I C HEAT RATIO * STAND- 

C ARD ACCELERATION OF GRAVITY • UNIVERSAL GAS CONSTANT). FOR EARTH, GASFAC= 

C 20.064881 (METERS / SEC**2 / KELVIN DEGREE t ** 1/2 . 

C 

COMMON C 
C 

DIMENSION A ( 600 ) » BI700), C(4000>, 

1VATM( 3 ) , P ( 3 ) , X IFT ( 3 ) , DRAG ( 3 ) , PAR(3>, XI 100) 

C 

EQUIVALENCE 


UA ,C 

l LIU, (ALPHA , A 

( 49)), (AREA 

B 

( 6) ) 

2 ( B , C 

(1111)1 , (BETA , A 

( 50)), (CD 

A 

( 165)) 

3 { CD I , A 

( 163)), (CL , A 

( 164) > , (COSALF 

B 

( 48) ) 

4 { CQSBET , B 

( 49)), (DNSITY,B 

( 29)), (DRAG 

B 

( 69) ) 

5 (GASFAC , A 

( 46)), (P , B 

( 84)), (PAR 

B 

( 60) ) 

6IPMAGN ,B 

( 50)), (Q , B 

( 59) ) , l R 

, B 

( 102)) 

7 ( S INALF , B 

( 46) ) , ( S INBET , B 

( 47)|,(TM 

, B 

( 34) ) 

8IVATM , B 

( 97)1, (VMACH , B 

( 38 ) ) , ( VQ 

r B 

( 100)) 

9 ( VQSQRD, B 

( 101)1, (X * B 

t 401)1, IXIFT 

, B 

( 72) ) 


C 

Q = 0.5*DNSITY*VQSQRD 
QVAL = Q*ARfcA/X(2) 

VMACH=SQRTF( VUSQRD/TM ) /GASFAC 
C 

C COMPUTE THE X,Y,Z COMPONENTS OF LIFT. 

IF (ALPHA) 2,1,2 

1 CL = 0.0 
C0I=0.0 
GO TO 4 

2 CL ■= QUAD! VMACH »2)*S1NALF 
AA = QVAL*CL/PMAGN 

AB = SINBET/VQ 
DO 3 K=1 , 3 

3 XIFT(K) = AA*(AB*PAR(K)+COSBET«P(K) ) 

7 CD 1= QUAD! VMACH,3)*CL*CL 

C 

C COMPUTE THE X,Y,Z COMPONENTS OF DRAG. 

4 CD = CDI+QUAD( VMACH, 4) 

AC = -CD*QVAL/VQ 

DO 5 K=l,3 

5 DRAG! K ) = AC*VATM(K) 

6 RETURN 
END 



FUNCTION ARCTAN <Y»X) 

C THE FORTRAN II LIBRARY ATANF ( + OR - Z*TAN ! THETA ) > USES A SINGLE 

C ARGUMENT WITH ITS SIGN TO GIVE THETA IN THE FIRST (+Z) OR FOURTH 

C (-Z) QUADRANT. 

C 

C THE ARCTAN FUNCTION MAY BE USED IF + OR - Z IS DERIVED FROM A 

C FRACTION SO THAT ARCTAN (Y,X) = TAN-1 ( ( +OR-Y=SlN( THETA) ) / ( +OR— X= 

C COS(THETA))). THUS THE ARCTAN (Y,X) GIVES THETA IN ITS PROPER 

C QUADRANT FROM -180 DEGREES TO +180 DEGREES. 

C 

IF (X) 2,1,2 

1 ARCTAN 3 S I GNF ( 1.57079632? Y) 

GO TO 4 

2 ARCTAN=ATANF( Y/X) 

IF ( X ) 3,1,4 

3 ARCTAN=ARCTAN+$IGNF( 3. 14159265, Y) 

4 RETURN 
END 


C 

C 

C 

C 

C 

c 

c 


c 


c 


SUBROUTINE C0NVT1 !V,AMC) 

THIS ROUTINE COMPUTES — (1) ANGULAR MOMENTUM, AMC(4) 

(2) ANGULAR MOMENTUM SQUARED, AMC l 5 > 

13) X, Y, Z COMPONENTS OF ANG. MOM., AMC(J) 

(4) VELOCITY, V ( 4 ) 

(5) VELOCITY SQUARED, V(5) 

COMMON C 


DIMENSION A ( 600 ) , BI700), C14000), 

1 AMC f 3 > , V t 5 ) , RB(3), IND(3> 

EQUIVALENCE 

1 1 A ,c ! IU),(B ,C Illll) ) , ( I NO , A ( 60) ), 

2 ( RB ,B ( 193)) 


1 


DO l J 1 = 1 , 3 
J2= IND I J 1 ) 

J3= I ND 1 J2 ) 

AMC { J3 ) = Rb(JL) *V ( J2 )— RB ! J2 ) * V ( J 1 ) 

AMC ( 5 ) = AMC t 1 ) **2+AMC ( 2 ) **2+AMC ( 3 ) **2 
AMC (4) 3 SQRTF { AMC ( 5 ) ) 

V( 5 ) = V!l)**2+V!2)**2+V!3)**2 
V ( 4 ) = SQRTF t V ( 5 ) ) 

RETURN 

END 


C 

c 

c 

c 

c 

c 

c 

c 


c 


c 


SUBROUTINE C0NVT2 

THIS ROUTINE CONVERTS RECTANGULAR COORDINATES INTO ORBIT ELEMENTS. 
RECTANGULAR COORDINATES- POSITION COMPONENTS , X, AND VELOCITY COMPONENTS , VX. 
THE ORBIT ELEMENTS ARE IN THE ORBELS ARRAY- 
ID ECCENTRICITY (4) INCLINATION 

(2) ARGUMENT OF PERICENTER 15) MEAN ANOMALY 

(3) LONGITUDE OF ASCENDING NODE 16) SEMILATUS RECTUM 

COMMON C 

DIMENSION A ( 600 ) , BI700), CI4000), 

1 AMC ( 3 ) , ORBELS [ 6 ) , RB(3) 


EQUIVALENCE 


1IA 

»C 

( 11) ) 

(AM 

,B 

! 90) ) 

, ( AMSQRD, B 

( 91) 

2! AMC 

, B 

I 87) ) 

( B 

,C 

( inn ) 

, (COSTRU.B 

( 53) 

3 ! EPAR 

, B 

( 26) ) 

( GK2M 

,B 

( 36) ) 

, (ORBELS, B 

( 116) 

4 ( R 

,B 

! 102)) 

, (RB 

, B 

( 193)) 

, ( S INTRU, B 

( 52) 

5 { TRU 

,B 

( 40) ) 

, ( V 

, B 

( 95) ) 

, ( VSQRD , B 

( 96) 

6! VX 

,B 

( 92) > 

, ( VY 

,B 

( 93)) 

,(VZ , B 

( 94) 


1 


2 

3 


4 

5 


6 

7 


ORBELS! 6) =AMSQRD/GK2M 

R=SQRTF(Rb! 1>**2+RB(2)**2+RB(3)**2) 

TRU=ARCTAN(AM/GK2M*(RB( 1)*VX*-RB!2)*VY+R8(3)*VZ) , ORBELS ( 6 ) -R ) 
IF(AMCIl)) 2,1,2 
ORBELS I 3 ) =0 ■ 

GO TO 3 

ORBELS ( 3 ) =ARCT AN I AMC 1 1 ) ,-AMC { 2 ) ) 

OR BELSI4)=ARC TAN (SQRTF! AMC ( 1 ) **2+AMC ( 2 ) * *2 ) , AMC ( 3) ) 
SN00E=SINF(URBELS(3)) 

CNODE=COSF(ORBELS(3) ) 

AA=R8( l)*CNUDE+RB(2)*SN0DE 

AB = RB( 3)*SINF(0RBElS(4) ) +COSF { ORBELS ( 4 ) ) * ( R0 ( 2 ) *CNQDE-RB I 1 ) *SNOD£ ) 
ORBELS! 2 ) =ARCTAN ! AB, AA ) -TRU 

ORBELS! 1 )=SURTF(ABSF( 1 . +ORBELS ! 6 ) • ! VSQRD/GK2M-2 . /R ) ) ) 

EPONE=SQRTFl l.+ORBELS! 1) ) 

E2Ml*l.— ORBELS! 1 ) **2 
EPAR=SQRTFIABSF(E2M1) ) 

SINTRU=SINF!TRU) 

COSTRU-COSF! TRU) 

EPAS-SQRTF ! ABSF ! 1.— ORBELS! 1) ) ) *S INTRU/ ( 1 . +COSTRU ) 

ETHETA=ORBELS! l)*SINTRU/( l.+ORBELS! 1 ) *CO$TRU ) *EPAR 
IF ( E2M 1 ) 5,6,6 

ORBELS! 5 ) = LUGF( ( EPONE+EPAS ) / ( EPONE-EPAS ) )-ETHE TA 
GO TO 7 

ORBELS! 5 ) *2. *ARCTAN! EPAS , EPONE ) -ETHETA 

RETURN 

END 
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c 

c 

c 

c 

c 

c 


c 


c 


c 

c 


c 

c 


c 

c 


SUBROUTINE ERRGRZ 

THIS SUBROUTINE COMPUTES THE RELATIVE ERRORS BETWEEN THE R-K AND LOW-ORDER 
INTEGRATION SLHEMES. IT ALSO COMPUTES THE ERROR COEFFICIENT, A, ANO SAVES 
THE ERROR DATA WHEN EREF HAS A - SIGN. THE BRANCH ON IMODE DETERMINES 
WHICH SET OF NORMALIZING FACTORS ARE TO BE USED. 

COMMON C 

DIMENSION A ( 600 1 * BC700), C(4000>, 

1 RELERR (8), XPRIM (200), XINC (100) 

EQUIVALENCE 

1 ( A ,C 

2 ( B ,C 

3 ( EREF , A 
4(R ,B 

5 ( V , B 

EQUIVALENCE 

E2 = 0. 

RELERRt 2 ) = X INC ( 2 ) /XPR IMI 2 ) 

IF < IMQDE-1) 2,1,2 


( 1 1 ) ) , ( A 1 , B ( 10)1, (A2 , B ( ID), 

( 1 1 1 1 ) ) , ( DELT , B < D),(E2 ,B ( 18)), 

{ 13)), (IMODE , A ( 1 ) ) , ( I NOERR, B ( 51)), 

( 102 ) ) t ( STEPGO, A ( 41 ) ) , ( STEPNO, A ( 42)), 

( 95)), (XINC ,8 ( 60D), (XPRIM ,C ( 711)) 

(RELERR, XINC) 


COMPUTE THE NORMALIZED INTEGRATION ERRORS FOR THE ORBIT ELEMENTS. 
1 RELERR! 3 )*X1NC( 3)/ (XPRIM(3)+1. ) / 10. 

RELERR! 8 ) -XINC ( 8 ) /XPR I M( 8) /10. 

DO 10 J= 1 , 4 

10 RELERRt J+3)=XINC( J+3)/ 62. 831853 
GO TO 3 


COMPUTE THE NORMALIZED INTEGRATION ERRORS IN RECTANGULAR VARIABLES. 

2 VI * V+100. 

DO 20 J= 1 , 3 

RELERR! J+2) =XINCt J+2 )/ VI 
20 RELERR ( J+5 ) =X INC ( J+5 ) /R 

SELECT MAXIMUM ERROR, COMPUTE ERROR COEFFICIENT, POSSIBLY SAVE ERROR DATA. 

3 DO 5 J-2,8 

IF { ABSF(RELERR( J) )-E2) 5,5,4 

4 K=J 

E2 = ABSF ( RELERR ( J ) ) 

5 CONTINUE 

E2 = E2 ♦ 2E-8 
A 1 = A2 

A2 = LOGF ( E2 ) — 5. *LOGF I ABSF ( DELT ) ) 

IF (EREF) 6,7,7 

6 WRITE TAPE 4 , STEPGO , STEPNO, XPR IM (l> , DELT , A2 , E2 ,( RELERR ( J ), J=2 , 8 >, K 
1NDERR = INDERR + 1 

7 RETURN 
END 


SUBROUTINE EQUATE 

C THIS SUBROUTINE IS CALLED FROM NBODY TO EVALUATE THE DERIVATIVES OF THE 

C VARIABLES OF INTEGRATION. EITHER RECTANGULAR COORDINATES OR ORBIT ELE- 

C MENTS MAY BE USED AS THE VARIABLES OF INTEGRATION, BUT IN THE CASE OF THE 

C LATTER, THE CORRESPONDING RECTANGULAR COORDINATES MUST FIRST BE FOUND. 

C THIS IS DONE AT THE BEGINNING THRU THE USE OF KEPLERS EQUATION. THE 

C PERTURBATING ACCELERATIONS ARE FOUND BY CALLING VARIOUS OTHER SUBROUTINES 

C ANO THEIR SUM RESOLVED ALONG THE X,Y,Z AXIS. FINALLY, THE DERIVATIVES 

C ARE CALCULATED. IN THE CASE OF ORBIT ELEMENTS, THE X,Y,Z PERTURBATING 

C ACCELERATION COMPONENTS MUST FIRST BE RESOLVED INTO C I RCUMFERENT I AL , RADI AL 

C AND NORMAL COMPONENTS. THIS ROUTINE ALSO CHANGES THE INTEGRATION VARI- 

C A8LES FROM URBIT ELEMENTS TO RECTANGULAR VARIABLES IF THE ECCENTRICITY 

C APPROACHES UNITY. THE X, XPRIM, AND XDOT ARRAYS ARE AS FOLLOWS. 

C 

C 

C X ORBIT ELEMENTS RECTANGULAR COORDINATES 

C 

C 1 

C 2 

C 3 

C 4 

C 5 

C 6 

C 7 

C 8 

C 

C 

C 

0 

1 
2 

3 

4 

5 XDOT (100) 

C 


TIME 

TIME 

MASS 

MASS 

ECCENTRICITY 

X-VELOCITY 

ARGUMENT OF PERICENTER 

Y-VELOCITY 

ARGUMENT OF ASC. NODE 

Z-VELOCITY 

INCLINATION 

X 

MEAN ANOMALY 

Y 

SEMILATUS RECTUM 

Z 


c 


IMENSION 

A! 600 ) , 

BI700), C ( 4000 ) , 



XPRIM 

1 100,2) , 

, VX (3), 

QX 

(3), 

RB 

(3), 

NEFMRS (8), 

X 

( 100) 

XPR I MB 

( 100,2) , 

, FORCE (3), 

XI FT 

(3) , 

DRAG 

(3), 

OBLAT (3), 

COMPA 

(3) , 
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equivalence 


1 ( A ,C 

( 111 ) , IAMSQRD.B 

( 91) ) 

, (ASYMPT, A 

( 7) ) 

2(B ,C 

( 1111) ) , (BNAME ,B 

I 122)) 

•(CINCL , B 

( 55) ) 

3 ( C I RCUM , B 

( 82 ) ) , ( COMPA ,B 

{ 63) ) 

, (CONSTU, A 

( 32) ) 

4 ( COSTRU, B 

( 53) } , I COSV , B 

( 57) ) 

, ( DRAG , B 

( 69) ) 

5 ( EMONE ,B 

< 2 8 ) ) , I EPAR ,B 

( 26) ) 

, ( E TOL , A 

( 30) ) 

6 ( E XMODE , B 

( 27)), (FLOW , B 

( 5) ) 

.(FORCE , B 

( 66) ) 

7 ( GK2M ,0 

( 36 ) ) , ( GKM ,B 

( 37) ) 

, ( IMUDE * A 

( 1) ) 

8 ( KSUB ,B 

( 19) ) , ( MBODYS, B 

( 42) ) 

, ( NEFMRS , B 

( 185)) 

9{NSTARF,B 

( 24 ) ) , ( OB L AT N , A 

( 40) ) 

, ( OBLAT ,6 

( 75) ) 

equivalence 

l ( PRESS ,B 

( 33 ) ) , ( P USHO ,B 

I 391)) 

, ( QX , B 

( 78) ) 

2 (RADIAL , B 

< 81) ) , IRATMUS.B 

( 23) ) 

, (RB , B 

( 193)) 

3 ( R ,8 

( 102) ) , I RSQRD ,B 

I 45) ) 

•(SINCL , B 

( 54) ) 

4( S INTRO, B 

( 52) ) , (SINV ,B 

( 56) ) 

,(SPU , A 

( 44)) 

5 ( TABLT ,8 

{ 20 ) ) , ( TOFF T , A 

( 24) ) 

, ( TRSEER, B 

( 8) ) 

6( TTEST , A 

( 54)), (U , A 

l 59) ) 

, t V , B 

( 95)) 

7 ( VSQRD ,B 

( 96)) , ( VX ,B 

( 92) 1 

, ( xoo r , b 

( 501)) 

8 ( X I F T ,B 

( 72 ) ) , ( XPR I M ,C 

( 711)) 

, (XPRIMB,C 

( 911)), 

9 ( X ,B 

( 401) ) , (ZORMAL,B 

( 83) ) 

, ( ZN , B 

( 43) > 


TAtiL T = X( ll/SPU+TOFFT 
IM0DE= IMODE 

1 GO TO 12,16, 16) , IMODE 
C 

C STATEMENTS 2 TO 16 FINO THE RECTANGULAR POSITION ANO VELOCITY FROM ORBIT 

C ELEMENTS ANO TRUE ANOMALY. THE TRUE ANOMALY IS FOUND FROM ITERATIVE 

C SOLUTION OF KEPLERS EQUATION. 

2 E2 = X( 3) **2 
E2M1 = l.-EZ 
EMONE = X ( 3 ) — 1 • 

EPAR=SQRTF(AB$F( E2M1 ) ) 

VC I RCL = GKM/ SQKTF ( X { 8 ) ) 

C 

C COMPUTE SINt ANO COSINE OF TRUE ANOMALY. 

C PART A. E= 1 

3 IF (EMONE) 10,4,5 

4 SINTRU * 0. 

COS TRU = 1. 

GO TO 14 

C 

C PART d. E IS GREATER THAN 1 

5 00 7 J=l»100 

DELM = X(7>-U+X(3)*SINHF(U) 

ECO SU = X( 3)»CDSHF(U) 

UELU = DELM/ ( l.O-ECOSU) 

U = U+DELU 

6 IF { ABSF ( DELM ) — CONSTU ) 9,9,7 

7 CONTINUE 
ASYMPT = 1.0 

IF (MdOOYS) 8,23,8 

8 CALL EPHMRS 
GO TO 23 

9 COSU = COSHF(U) 

DEMI = 1 .-X l 3) *CQSU 
COSTRU = ( CuSU-X ( 3 ) ) / DEM 1 
SINTRU =-EPAR*SINHF( U) /DEMI 
GO TO 14 

C 

C PART C. E IS LESS THAN 1 

10 DO 12 J=l, 15 

DELM = X(7)-U+X(3)*SINF(U) 

ECOSU = X(3)*C0SFfU) 

DELU * DELM/ ( l . O-ECGSU+O . 0 1 *ECOSU* * 3 ) 

U = U+DELU 

11 IF l AtJSF(UELM)-CONSTU) 13,13,12 

12 CONTINUE 

WRITE OUTPUT TAPE 6,55, U, DELU 
CALL EXIT 

13 COSU = COSF(U) 

0EM1 * l.-Xl 3 ) *COSU 
CUSTRU = ( CuSU-X ( 3 ) J/DEM1 
SINTRU = EPAR»S1NF(U) /DEMI 

14 PDVR = l.+XI 3)*C0STKU 
C 

C COMPUTE POSITION AND VELOCITY FROM ORBIT ELEMENTS AND TRUE ANOMALY. 

C ALSO, CLEAR THE PERTURdAT I NO ACCELERATIONS. 

15 SOMEGA = S INF( X( 4) ) 

COMEGA = COiF ( X ( 4 ) ) 

SNOOE = S INF ( X ( 5 ) ) 

CNCUE = COSF ( X ( 5 ) ) 

SINCL = S I NF I X ( 6 ) ) 

CINCL = COSM X ( 6 ) ) 

S I NV=SINTRU*CUMEGA+COSTRU* SOMEGA 

COS V = COS TRU » COMEGA- S I NTRU* SOMEGA 

AR=COSV*CNUuE-SINV*SNUDE*CINCL 

81=SINV*CN0uE+C0SV*SN0DE*CINCL 

C1=COSV*SNODE+SINV*CNODE*C INCL 

D1=S INV*SNOUE— COSV+CNODE+C INCL 

El = X( 3>*S0MEGA+SINV 

FI = X( 3>*CUMEGA+C0SV 

AS = £ 1 *CN0UE + F l » SNOOE *C INCL 

02 = F 1 »CNDUE *C I NCL— E 1 * SNOOE 

R = X ( 8 ) /PD VR 

RSQRD = R*R 

SINVY = S1NV*S INCL 

RB ( 1 ) = K* AR 

RB ( 2 ) = R«Cl 

RB ( 3 ) = K*S INVY 

VX( l ) VCIRCL»AS 

VX(2)=VCIRCL*B2 

VX(3)=VCIRCL*F1*SINCL 

GO TO 18 
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c 


16 DO 17 K=1 i 3 

VX ( K ) = X(K+2) 

17 RB(K> = X ( K + 5 ) 

RSQRD = RB ( 1 ) *RB ( 1 ) + RB(2)*RB(2) + RB(3)*RB(3) 

R= SQRTF (RSQRD ) 

18 VSQRD=VXtl)#VX( l)+VX(2)*VX(2)+VX(3)»VXI3) 

V = SQRTF1 VSQRD) 

DO 19 1=1,15 

19 FORCE( I) = 0. 

C 

C TEST FOR PRESENCE OF PERTURBING BODIES* 

IF (MBQDYSl 20,21,20 

20 CALL EPHMRS 

21 IF (XABSFI IMOUE)-l) 26,22,26 
C 

C TEST FOR CHANGE FROM ORBIT ELEMENTS TO TEMPORARY RECTANGULAR 

C COORDINATES IF E IS TOO NEAR TO UNITY. 

22 IF ( ETOL— ABSF ( EHONE ) ) 26,23,23 

23 IF (IMODE) 54 , 24,24 

24 IMODE=— 3 

IF INSTART) 25,54,25 

25 TTEST = X(l) 

27 CALL TESTTR 

C 

C TEST FOR OBLATENESS PERTURBATION COMPUTATION. 

26 IF I OBLATN— BNAME >30,29, 30 

29 CALL OBLATE 
C 

C TEST FOR PRESENCE OF THRUST. 

30 XD0TI2) = -FLOW 

IF (R-RATMOS) 31,31,32 

31 CALL ICAO 
GO TO 33 

32 PRESS=0. 

33 IF (PUSHO) 37,36,37 

36 ASSIGN 40 TO NDONE 
GO TO 38 

37 CALL THRUST 
ASSIGN 41 TU NOONE 

C 

C TEST FOR EXISTENCE OF ATMOSPHERE- FIND AERODYNAMIC FORCES. 

38 IF IPRESS ) 39,42,39 

39 GO TO NDONE, (40,41) 

40 CALL THRUST 

41 CALL AERO 
C 

C SUM COMPONENTS OF THE PERTURBING ACCELERATION. 

42 DO 43 J= 1 , 3 

43 COMPAIJ) = -QX ( Jl+OBLATIJ) + FORCE (J)+XIFTtJ)+DRAG(J) 

44 GO TO 147,46,45) , IMODE 
C 

C COMPUTE DERIVATIVES FOR THE RECTANGULAR VARIABLES OF INTEGRATION. 

45 AA = GK2M/R/RSQRD 
DO 46 K= 1 , 3 
XD0TIK+5) = XlK+2) 

46 XDOTtK+2) = COMPA ( K ) -AA*X1 K+5 ) 

GO TO 54 

C 

C COMPUTE THE DERIVATIVES OF THE ORBIT ELEMENTS. (AFTER RESOLVING 

C PERTURBATINb ACCELERATION INTO CIRCUMFERENTIAL, RADIAL, NORMAL COMPONENTS) 

47 CIRCUM = COMPA( 3 ) *C0S V*S I NCL-COMPA { l ) *B 1 -COMPA ( 2 ) * D1 
RAD I AL = COMPA ( l ) » AR+COMPA { 2 ) *C 1+COMPA ( 3 ) *S I NV Y 

ZORMAL = COMPA { 1)«SN0DE*SINCL— C0MPA(2) *CNODE*S I NCL+COMPA ( 3 ) *C INCL 
ZN=VCIRCL*E2Ml*EPAR/X( 8) 

RDVPP1 = l./PDVR + 1. 

RDVA = E2M1/PDVR 

XDOT ( 8 ) = 2.*R/VCIRCL*C IRCUM 

IF I X ( 3 ) ) 48,48,49 

48 CSQRD = C IRCUM*C I RCUM 
RASQRD = RAuIAL-RADIAL 

DEMI = (4.*LSQRD+RASQRD)*VCIRCL 
C 

C TEST FOR IN-PLANE PERTURBATION. 

IF (DEMI) 57,56,57 

56 X00TI3) = 0. 

XDOT ( 4 ) = 0. 

XDOT ( 7 ) = 0. 

GO TO 50 

57 V0V2R=VC IRCL/R/2 . 

XDOT ( 3 ) = SQRTF ( 4. «CSQRD+RASQRD ) / VC I RCL 
XDOT 14) = VUV2R+(2.*CSQRD+RASQRD)/DEM1*RADIAL 
XDOT ( 7 ) = ZN-VDV2R+(6.*CSQRD+RASQRD)/DEM1*RAD1AL 
GO TO 50 

49 XDOT ( 3 ) = (SINTRU*RADIAL+(PDVR-RDVA)/X(3)*CIRCUM)/VCIRCL 
XD0T(4) = (SINTRU/X(3)*RDVPP1*CIRCUM-C0STRU*RA0IAL/X(3) )/ VC I RCL 
XDOT ( 7 ) =ZN+EPAR/ VC I RCL * ( 1 COSTRU/ X ( 3 ) —2 ./PDVR)*RADIAL“(S1NTRU/X(3)* 

IRDVPP 1*C IRCUM) ) 

50 IF(SINCL) 51,52,51 

51 XDOT ( 5 ) = SINV/ S INCL*ZORMAL/ VC IRCL/PDVR 
GO TO 53 

52 XDOT ( 5 ) = 0. 

53 XDOT ( 6 ) = CUSV* ZORMAL/ PD VR/ VC I RCL 

54 RETURN 

55 F0RMAT(41H0KEPLERS EQUATION CONVERGENCE FAILURE, U=G15.8,7H DELU= 

IG15.8) 

END 
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subroutine ephmrs 

C SUBROUTINE EPHMRS IS CALLED TD COMPUTE THE POSITIONS OF THE PERTURBING 
C BODIES RELATIVE TO THE VEHICLE AND, FROM THESE, THEIR PERTURBING ACCELERA- 
C TIONS UPON THE VEHICLE. OCCASIONALLY THIS ROUTINE IS CALLED FOR THE PURPOSE 
C OF TRANSLATING THE ORIGIN IN WHICH CASE (TRSFER=1) THE RELATIVE VELOCITIES 
C ARE ALSO CALCULATED. IF A BODYS POSITION IS TO BE COMPUTED FROM AN ELLIPTIC 
C APPROXIMATION SUBROUTINE ELIPSE IS CALLED. OTHERWISE, THE POSITION WILL BE 
C CALCULATED IN EPHMRS FROM THE PRECISION TAPE EPHEMER IS. THE DO 19 LOOP 

C ENCOMPASSES ALMOST THE ENTIRE EPHMRS SUBROUTINE AND ,IN EFFECT, ELIPSE TOO. 

C 

COMMON C 


C 


c 


c 

c 


c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


c 

c 


c 

c 

c 

c 


c 

c 


c 

c 

c 


c 

c 

c 


DIMENSION A1600I, 8(700), CI4000), 

1 QX( 3) , IbODYtB), EFHRS ( 7 ) , XP(3,8), RB(3,8), R (8), TIM(7), 

2 NEFMRS 1 8 ) , TDAT A( 6,3,7) , TDEL ( 7 ) , BHASS<8), VEFM(3,8), DATA ( 21 ) 

3 , TOATI 18,7) 


EQUIVALENCE 
11A ,C ( 

2IBMASS , B 1 

3(1 BODY ,B ( 

4 ( QX , B ( 

5 ( SQRDK , B ( 

6 ( TOATA ,0 l 

7 ( TRSFER, B ( 

EQUIVALENCE ( IBF 


1 1) ) , ( AU , A 

137) ) , ( DTOFF J , A 
177 ) ) , ( MtJODYS, B 
78)), (R ,B 

35) » , (SPD , A 
265)1, (TDEL ,B 
8) ) , ( VEFM ,B 
, F IB ) , ( TDAT , TDATA ) 


( 29)), (B ,C 
( 23) 1 , (EFMRS ,B 
( 42) ), (NEFMRS, B 
( 102)), (RB ,B 
( 44) ), (TABLT ,B 
( 170)), (TIM ,B 
( 241)), (XP , B 


(HID) 
( 130)) 
( 185)) 
( 1931) 
( 20 ) ) 
( 163)) 
( 217)) 


PART 2. SET INDEXS, FIND POSITION IF ELLIPSE IS USED (NEFMRS = 20 OR UP). 

DO 19 JB=1,MB0DYS 

JB1 = JB + 1 

IBF = IBODYl JB1 ) 

IB = XABSFI IBF) 

IF (NEFMRS! JB)-20) 2,2,1 
1 CALL ELIPSE (JBL) 

IF (TRSFER) 12,12,17 


PART 3. TAPE EPHEMER I S IS TO BE USED. FIND DIFFERENCE IDT) BETWEEN 
CURRENT PROBLEM TIME ( DT0FFJ+TA8LT > AND MIDPOINT TIME (TIM) OF CURRENTLY 
STORED TAPE DATA. THEN SEE IF CURRENT DATA IS OKAY. TDEL = TIME INTERVAL 
ON EITHER SIDE OF TIM FOR WHICH CURRENT DATA IS GOOD. 

2 DT = TABL T - (TIM(JB) — DTOFF J ) 

IF !ABSF!DT)— TDEL(J8) ) 10,10,3 


PART 4A . CURRENT DATA NOT OKAY. READ IN NEXT DATA SET. IF DT IS 
BACK UP THE TAPE 2 RECORDS BEFORE READING. 

3 IF (DT) 4,5,5 

4 BACKSPACE 3 
BACKSPACE 3 

5 READ TAPE 3, ( DATA ( I ) , 1=1,21) 

PART 4B. IF THIS DATA IS FOR A BODY IN THE BNAME LIST, STORE IT. 

(IF NOT STORED, WE MIGHT HAVE TO RETURN FOR IT.) IF ELLIPSE OATA IS 
PROVIDED FOR THE BODY FOUND, BY-PASS THE TAPE DATA AND READ IN NEXT SET. 

DO 7 J = l , MBODYS 

IF ( ( DATA! D+EFMRSIJ) )*(— (DATAll) * EFMRS ( J ) ) ) ) 7,6,7 

6 IF ( NEFMRS ( J )— 20 ) 8,8,3 

7 CONTINUE 
GO TO 3 

PART 4C. MOVE THE DATA INTO PLACE AND THEN GO BACK AND SEE IF IT IS OKAY. 

8 TIM(J) = DATA ( 2 ) 

TDEL ( J ) = DATA ( 3 ) 

DO 9 J J = 1 , 18 

TDAT ( J J , J ) = DATA(JJ+3) 

.9 CONTINUE 
GO TO 2 


PART 5. CURRENT DATA IS OKAY. GET POSITION FROM THE POLONOMIAL 
P = A + BX + CX**2 + DX**3 + EX**4 +■ FX**5. 

10 DO 11 K= 1 , 3 

XP(K.JBl) = TDATA { 1 , K , JB ) 

DO 11 K T = 2 , 6 

XP ( K , J B 1 ) = XP ( K , JB1 ) * DT +TDATA ( KT , K, J B 1 

11 CONTINUE 

IF (TRSFER) 12,12,15 

PART 6. COMPUTE DISTANCE FROM REFERENCE AND FROM ROCKET . 

12 DO 13 K=I,3 

XP1K.JB1) = XP ( K , I B ) +XP(K,JB1)*SIGNF(AU,FIB) 

13 RB ( K , JB1 ) = RB ( K , 1 ) - XP(K,JB1> 

PART 7. COMPUTE PERTURBING ACCELERATIONS (QX). 4194304=2**22 IS REMOVED 

TO PREVENT OVERFLOW. 2048=2**11 AND 8589934592=2**33 RESTORE THE SCALE. 
PRSQRO = ( RB ( 1 , JB1 ) **2 ♦ RB(2,JB1)*»2 + RB( 3, JB1 ) **2 ) /4194304. 

RRELL = SQRTF(PRSQRD) 

RSQRD = ( XP ( 1 , JB1 ) **2 + XPI2,JB1)**2 + XP ( 3, JB1 ) **2 ) /4194304. 

RCUBE = RSQRD * SQRTF ( R SQRO ) 

PRCUBE = PRSQRD * RRELL 
R(JB1) = RRfcLL*2048- 
DO 14 K= 1 , 3 

14 QX ( K ) =SQRDK • BMASS(JBl) * (( XP ( K , JB1 ) /RCUBE ) + RBI K , JB 1 ) /PRCUBE ) / 

1 8589934592. + QX(K) 

GO TO 19 

PART 8. COMPUTE VELOCITY FROM V = B + 2CX + 3DX**2 + 4EX**3 + 5FX**4 
AND FROM REFERENCE BOOY VELOCITY ( VEFM ( IB)). 

15 DO 16 K= 1 , 3 
VEFM(K.JBl) = 0. 

DO 16 KT = 1 » 5 

16 VEFM ( K, JB1 ) = ( VEFMtK, JB1) *DT+TDATA ( KT , K, JB ) *FLOATFl -KT+6 ) ) 

17 DO 18 K=1 , 3 

18 VEFM( K, JB1 ) = VEFM( K, I B ) + V£FM(K»JB1)*S IGNF (A U/ SPD, FIB) 

GO TO 12 

19 CONTINUE 
RETURN 
END 
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SUBROUTINE EXTRA 
C 

C THIS ROUTINE IS EXECUTED BETWEEN FLIGHTS AND MAY THEREFORE BE EXPANOED TO 

C DO ADDITIONAL COMPUTATION BETWEEN SUCCESSIVE FLIGHTS. 

C 

COMMON C 
C 

DIMENSION A1600), B1700), CI4000) 

C 

EQUIVALENCE 

1(A ,C ( 11)1, (B ,C (1111)), (QMAX , B ( 44)), 

2 ( S I GNAL * B ( 31)) 

C 

SIGNAL * 0. 

QMAX - 0. 

RETURN 

END 


SUBROUTINE EXTRAS 

C THIS ROUTINE IS EXECUTED BETWEEN STAGES AND MAY THEREFORE BE EXPANDED TO 

C DO CALCULATIONS BETWEEN SUCCESSIVE STAGES OF A FLIGHT. 

C 

RETURN 

ENO 


SUBROUTINE ELIPSE (JBI) 

C THIS SUBROUTINE IS CALLED FROM EPHHRS TO COMPUTE THE POSITION OF A BODY 

C USING APPROXIMATE ELLIPTIC DATA. THE VELOCITY IS ALSO COMPUTED IF THE 

C ORIGIN IS BEING TRANSLATED I TRSFER= I .0 ) • THE ELLIPSE DATA IS READ FROM 

C INPUT CARDS AND ORGANIZED IN SUBROUTINE ORDER. TPD IS TIME SINCE PERIHELION 

C PASSAGE, ZM IS MEAN ANOMALY, U IS ECCENTRIC ANOMALY. 

C TDATA ARRAY - ( K ) SEMILATUS RECTUM (K+7) PERIOD 

C t K+l ) ECCENTRICITY (K+8) SIN OMEGA 

C ( K+2 ) OMEGA (K + 9) SIN NODE 

C ( K+3) NODE (K+10) SIN INCLINATION 

C (K+4) INCLINATION (K+ll) COS OMEGA 

C (K+5) JD OF PERIHELION (K+12) COS NODE 

C I K+6 ) FRACTIONAL PART OF (K+5) (K+L3) COS INCLINATION 

C 

COMMON C 
C 

DIMENSION A ( 600 ) , B(700), CI4000), 

I XP ( 3 , 8 ) , VEFM ( 3,8), TUATAl 121 ) 

C 

EQUIVALENCE 

1 ( A ,C ( ID), IB ,C ( 1111) ) , (CONSU ,A ( 311), 

2 ( DTOFF J , A ( 23 ) ) , ( TABLT ,B ( 20)), (TDATA ,B ( 265)), 

3 ( TRSFER, B ( 8)), (VEFM ,B ( 24l)),(XP ,B ( 217)) 

C 

K =■ 18* ( JBI— 2) + l 

TPD = ( DTOFF J-TDATA( K+5) )+( TABLT— TDATA ( K+6 ) ) 

ZN = 6. 28318533/ TDATA ( K+7 ) 

ZM * ZN+MODF (TPD, TDATA! K+7) ) 

C 

C GET THE SINE(SINTRU) AND THE COSINE (COSTRU) OF THE TRUE ANOMALY 

C BY ITERATING KEPLERS EQUATION. THEN COMPUTE X,Y,Z (XP). 

U = ZM+TDATA(K+1)*SINF(ZH)+.5*TDATA(K+1)**2*SINF(2.*ZM) 

DO 1 J=i,10 

DELM = ZM— U+TDATA( K+1)*SINF(U) 

DELU = DELM/ ( 1 .-TDATA I K+1)*C0SF(U) ) 

U - U+DELU 

IF ( ABSF(DELM)-CONSU) 2,2,1 

1 CONTINUE 

2 COSU = COSF(U) 

DENOM =* l.-TOATA(K+l)*COSU 
COSTRU = ( COSU— TDATA ( K+l ) ) /DENOM 
R - TDATA ( K) / ( 1. + TDATA (K+l ) * COSTRU) 

SINTRU * SQRTFl 1.— TDATA (K+l) **2 )*S INF (U) /DENOM 
SINV * SINTRU*TDATA( K+ll )+CQSTRU*TDATA( K+8) 

COSV = C0STRU*TDATA(K+ll)-SlNTRU*TDATA(K+8) 

XPIltJBl) = R*(CDSV*TDATA(K+12 ) — S INV+TDATA ( K+9 ) • TDATA t K+13 ) ) 

XP ( 2, JBI ) - R»(COSV*TOATA( K+9 ) +S INV*TDATA ( K+ 12 )* TDATA ( K+13) ) 

XP ( 3 , JB 1 ) = R*SINV*TDATA(K+10) 

IF (TRSFER) 3,4,3 
C 

C COMPUTE THE VELOCITIES FOR THE TRSFER OF ORIGIN. 

3 EX * TDATA(K+l)*TDATA(K+8)+SINV 
FX » TDATA(K+1) *TDATA( K+ll )+COSV 

CFACT = ZN*TDATA(K)/ (SQRTFt ( 1.— TDATA(K+1)**2)**3) ) 

AX * EX*TDATA(K+12)*FX*TDATA(K+9) * TDATA! K+13) 

BX - FX*TDATA(K+12)*TDATA( K+13) — EX*TDATA( K+9 ) 

VEFM ( 1 , JB 1 ) - — AX*CFACT 
VEFHI 2, JBI ) * BX+CFACT 
VEFM(3,JB1) * FX*CFACT*TDATA l K+10 ) 

4 RETURN 
END 
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REM SUBROUTINE EXADD ( A , B , C ) 

REM THIS ROUTINE WILL ADO IN DOUBLE PRECISION A QUANTITY C TO THE DOUBLE 
REM PRECISION VARIABLE A+B WHERE A IS THE MOST SIGNIFICANT PART AND B IS 
REM THE LEAST SIGNIFICANT PART. 



ENTRY 

EXADD 


COMMON 

-206 

Q 1 

COMMON 

1 

Q2 

COMMON 

1 

TEMPI 

COMMON 

l 

TEMP2 

COMMON 

HTR 

1 


BCI 

1, EXADD 


SXD 

—4,1 


SXD 

*— 4 , 2 

EXADD 

SXD 

—4,4 


CLA* 

1.4 


FAD* 

3,4 


STQ 

Q1 


FAD* 

2,4 


STQ 

Q2 


FAD 

Q1 


STQ 

Qi 


STO 

TEMPI 


CLA 

Ql 


FAD 

Q2 


STO 

TEMP2 


FAD 

TEMP2 


FAD 

TEMPI 


STQ 

Ql 


FSB 

TEMP2 


STO* 

1.4 


STQ 

Q2 


CLA 

Ql 


FAD 

Q2 


STO* 

2,4 


TRA 

END 

4,4 


SUBROUTINE ICAO 

C SUBROUTINE ICAO DETERMINES THE ATMOSPHERIC TEMPERATURE, PRESSURE, AND 
C DENSITY AS A FUNCTION OF ALTITUDE ABOVE THE EARTH IN ACCORDANCE WITH 

C THE 1962 U.S. STANDARD ATMOSPHERE (ICAO TO 20 KM.). A SHORT FAP 

C PROGRAM FOLLOWS ICAO WHICH PROVIDES A MEANS OF LOADING DATA INTO MACHINE. 

C IT MUST BE LOADED DIRECTLY AFTER ICAO. IF THE LENGTH OF ICAO IS CHANGED, 

C THE DATA MUST BE RELOCATED. 

C 

C R IS 

C ALT IS 

C TABLE H IS 

C 

C ALM IS 

C TMR IS 

C REF P IS 

C TM IS 

C 

C PRESS IS 

C DNS I T Y IS 

C HEIGHT IS 

C 

COMMON C 
C 

DIMENSION A ( 600 ) , B1700), C(4000), 

I TABLEHU3), TMR123), REFP(23), ALM(23) , RB(3> 

C 


EQUIVALENCE 


1 ( A ,C 

( 

11) ) , (ALT 

, A 

( 4 ) ) , ( B ,C 

( 1111) ) 

2 ( DNS I T Y , B 

( 

29) ) , 


(OBLATN, A 

( 40) ) 

3 ( PRESS t B 

( 

33) ) , (R 

, B 

( 102)), IRB , 8 

( 193)) 

4 ( RE , A 

( 

25) ) , I TABLT 

, B 

( 20)), (TM , B 

( 34) ) 

5 ( RESQRD, B 

( 

7) ) 




EQUIVALENCE 

(TA0LEH(24), TMR) 

( TABLEHI 47 ) , ALM } , ( TABLEH{ 70 ) 

, REF P ) 


C 

IF ( OBLATN ) 102,101,102 

101 ALT = R - RE 
GO TO 103 

102 ALT = R— 6356783- 2B/SQRTF ( . 9933065783+. 00669342 16851 RB(3)/R)**2) 

103 IF (ALT-90000.) 105,104,104 

104 HEIGHT = ALT 
GO TO 106 

105 HEIGHT = ALT/ ( 1 • O+ALT/ 6356766. ) 

106 K=K 
C 

C FIND THE HEIGHT IN A TABLE OF BASE DATA. DATA ARE 

C ARRANGED IN DECENDING ALT WITH 21 REGIONS. ABOVE THAT, PRESSURE AND 

C DENSITY ARE SET = 0. TEMPERATURE IS SET TO 3000. 

1 IF ( K— 22 1 2,6,6 

2 IF l HE IGHT— TABLEHI K+l) ) 5,3,3 

3 K = K+l 
GO TO 1 

4 K - K— 1 

5 IF (K) 7,7,6 

6 HINC = HEIGHT - TABLEH(K) 

IF (H INC) 4,8,8 

7 K - 1 

8 IF (ALMIKM 9,100,9 
C 


DISTANCE TO CENTER OF EARTH IN METERS. 

VEHICLE ALTITUDE ABOVE EARTH IN METERS. 

METERS OF ALTITUDE FROM THE EARTHS SURFACE AND IS 
THE ARGUMENT OF ATMOSPHERE PROPERTY TABLE. 

THE MEAN SLOPE OF THE TABLE H VS. TM CURVE AT TABLE H. 

TM AT TABLE H. 

THE PRESSURE IN MILLIBARS AT TABLE H. 

THE TEMPERATURE TIMES STD. MOLECULAR WEIGHT / ACTUAL 
MOLECULAR WEIGHT. DEGREES KELVIN. 

PRESSURE IN MILLIBARS. 

DENSITY IN KILOGRAMS PER CUBIC METER. 

EITHER GEOPOTENT I AL ALTITUDE OR GEOMETRIC ALTITUDE IN METERS. 
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c CONTROL COMES HERfc FOR NON I SOTHERMAL LAYERS 

9 TM * T MR { K ) + ALM(K)*H INC 
IF ( AL T— 90000 * > 107, 107,108 

107 PRESS * REFP ( K ) * ( T HR (K)/TM)*«{ . 034 163 1947/ ALM { K ) ) 

GO TO 10 

108 IF (K-KC) 109,110,109 

109 KC * K 

Cl * RE+TA8LEH(K> 

C2 * TMR < K ) / ALM ( K) 

C3 « l./<Ci-C2) 

C4 = -.0341631947*RESQRD*C3/ALM(K) 

110 PRESS » REFP(K) »EXPF ( C4* ( C 3* LOGF { C 1* ( H INC/C2+1 . ) / < RE+HE IGHT ))- 

1 HINC/C l/IRE+HE IGHT ) ) ) 

10 ONSITY » PRfcSS/TM/2. 87053072 
GO TO 13 

C 

C CONTROL COMES HERE FOR ISOTHERMAL LAYERS 

100 IF IK-22) 11,12,12 

11 TM = TMRIK) 

PRESS « R£FP(KI*EXPF(~.0341631947*HINC/TMR(K> ) 

GO TO 10 
C 

C CONTROL COMES HERE FOR EXTREME ALTITUDES 

12 PRESS = 0.0 
DNS I TY » 0.0 
TM = 3000. 

13 RETURN 
END 


REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

ORG 

Ai DEC 

DEC 
DEC 
DEC 

A2 DEC 

DEC 

A3 DEC 

DEC 
OEC 

A* DEC 

DEC 
DEC 
END 


THIS IS THE FAP PROGRAM WHICH LOADS ICAO DATA INTO MACHINE. 

THE 256 IN ORG 256 WAS FOUND BY SUBTRACTING 22 FROM THE OEC LOCATION 
OF REF P (FROM FAP LISTING OF ICAO, THIS WAS FOUND TO BE 278). 

THUS, 278-22*256. DISCARD THE FIRST TWO BINARY CAROS AFTER ASSEMBLY 
AND PLACE REMAINING CARDS IMMEDIATELY BEHIND ICAO BINARY DECK. 

Al IS REF P ( 23) 

A2 IS ALM ( 23) 

A3 IS T MR ( 23) 

A4 IS TABLE H(23) 

256 

0., 1. 19l8E-9,3.4502E-9, l . 0957E-8, 4. 0304E-8 , 1.0838E-7 
6.9604fc-7,1.6852£-6,2.7926E-6,3.6943E-6,5.06l7E-6, 2.5217E-5 
7. 3544t~5, 3.0075E-4, 1.6438E-3, .010377, . 182099,. 590005 
1. 10905,8.68014,54.7487,226.320, 1013.25 

0.»0. , .00 11,. 00 17,. 0026, .0033,. 004, .005, .007, .01,. 01 5, .02,. 01 

.005, .003,0., -.004, -.002,0. ,.0028, .001 ,0.,-. 0065 

0., 2700. 65, 2590. 65, 2420. 65, 2 160. 65, 1830 .65, 1 550. 65 , 1350. 65 

1210.65, 1110.65,960.65, 360.65,260.65,210.65, 180.65,180.65 

252.65,270.65,270.65,228.65,216.65, 216.65,288.15 

1E30, 7t5,6E5,5E5,4E5,3E5,2.3E5, 1.9E5, 1.7E5, 1 . 6E5, 1 . 5E5, 1 . 2E5 

1.1E5, IE5,.9E5, 79000. ,61000. ,52000. ,47000., 32000. ,20000. 

1 1000 . , 0 . 

1 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE NBODY 


NBOOY COMPUTES THE TRAJECTORY IN EITHER ORBIT ELEMENTS OR RECTANGULAR 
COORDINATES USING THE RUNGE-KUTTA TECHNIQUE. A LOWER ORDER INTEGRATION 
TECHNIQUE IS ALSO PERFORMED TO FACILITATE AUTOMATIC STEP SIZE CONTROL. 
THE X, XPRIM , XOOT , X I NC , ETC. ARRAYS ARE AS FOLLOWS. 


X 


ORBIT ELEMENTS 


RECTANGULAR COORDINATES 


TIME 

MASS 

ECCENTRICITY 
ARGUMENT OF PERICENTER 
ARGUMENT OF ASC. NODE 
INCLINATION 
MEAN ANOMALY 
SEMILATUS RECTUM 


TIME 

MASS 

X-VELOCITY 

Y-V6L0CITY 

Z-VELOCITY 

X 

Y 

Z 


IMOOt VARIABLES 

1 ORBIT ELEMENTS 

2 RECTANGULAR 

3 RECTANGULAR TEMPORARY 

4 EARTH SPHERICAL— CHANGE TO RECTANGULAR 

-1 ORBIT ELEMENTS — CHANGE TO RECTANGULAR 
-2 RECTANGULAR — CHANGE TO ORBIT ELEMENTS 

-3 ORBIT ELEMENTS — CHANGE TO TEMPORARY RECTANGULAR 
-4 EARTH SPHERICAL — CHANGE TO ORBIT ELEMENTS 


C 


c 


COMMON C 


DIMENSION 

1 XPRIM 

2 X 

3 XDOT 

4 AMC 

5 X WHOLE 


A (600) , 
( 100 , 2 ) 
( 100 ) , 

( 100 ), 

< 3) , 

( 6 ) , 


B( 700) , C ( 4000 ) , 

XPRIMB (100,2), 
XINC (100), 

RB (3), 

AK (3) , 

VX (3), 


XOOTPM (100,2), 
QLDINC (100), 

XK (100), 

AW (4), 

BEX (14) 
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EQU I VALENCE 


1 1 A 

,C 

I lin.iAi , b 

I 10) ) , I A2 , B 

( 11) ) 

2IAC0EF1 

,8 

I 12 ) ) , I AC0EF2 , 8 

I 13) ) , ( ACDEF3, B 

[ 14) ) 

3 1 AK 

, A 

I 51 ) ) , I AMC ,B 

( 87 ) ) , { AMSQRD, B 

( 91) ) 

41 AM 

,8 

( 90) ) , ( ASYMPT , A 

I 7 ) ) , 1 AW , A 

( 55) ) 

5 1 B 

, c 

(1111) ), I CONSTU , A 

I 32 ) ) , ( DEL T ,B 

( 1) ) 

6IOONE 

,B 

I 39)), (E2 , B 

I 18) ) , ( EMONE , B 

( 26) ) 

7 4 ERL I M T 

, A 

l 14) ), IETOL , A 

I 30) ), (EXMODE, 8 

I 27)) 

81GK2M 

, B 

I 36)), (H2 , B 

I 15)), (IMODE , A 

( 1) ) 

91 INOERR 

, B 

{ 51)), (KSUB , B 

I 19) ) , I MBODYS, B 

( 42) ) 

EQUIVALENCE 




1 1 NEQ 

, A 

( 2 ) ) , I NST ART , B 

I 24) ) , I OLDOEL, B 

( 9) ) 

2IQMAX 

, B 

I 44)), (RATIO , B 

( 58)), (RB ,B 

( 193)) 

3IREVS 

, A 

( 48)), IR , B 

I 102) ), (STEPMX,A 

( 16) ) 

4ISTEPG0 

.A 

I 41 ) ) , ( STE PNO , A 

I 42) ), (TRSFER, 8 

( 8)) 

5ITRU 

»B 

I 40) ) , I TTEST , A 

I 54 ) ) , ( VSQRD ,B 

( 96) ) 

6 1 VX 

, B 

I 92 ) ) , I XOOT ,B 

1 50 1 ) ) , ( X INC , B 

( 601)) 

7IXPRIM 

,c 

{ 711)),! XPR I MB, C 

( 911) ), (XWHOLE, B 

( 110)) 

8 I X 


I 40 1 ) ) , I ERLOG , B 

( 17) ) , ( EREF , A 

( 13) ) 

9 I Q 

, B 

I 59) ), (OUTPOT, B 

( 399)) 



C 

C PART 1. SET UP THE STARTING SEQUENCE FOR ERROR CONTROL AND DELAY CHECKING 

C THE ERROR UNTIL TWO STEPS ARE COMPLETED. THE ASSIGNED GO TOS NSTART AND 

C IBEGIN CONTROL STARTING. 

NEQ = NEQ 

1 DO 2 J=l,NEg 

XPR I M { J t 2 ) = XPRIM(J,I) 

XPR I MB 1 J , 2 J = XPRIMB(J*I) 

2 XIJ) = XPR IM I J » 1 ) 

NSTART = 0 
TRSFER = 0. 

H2 - DELT 
DELT = DELT/2. 

220 CALL EQUATE 

IF (OUTPOT) 222,221,222 

221 CALL OUTPUT 

222 DO 3 J= 1 , 3 
XWHOLE ( J I = VX( J ) 

3 XWHOLE ( J+3) = RB ( J ) 

C CHANGE INTEGRATION VARIABLES IF IMODE IS -. 

IF [IMODE) 4.5,5 

4 CALL TESTTR 
GO TO 1 

5 CALL TESTTR 

IF (TRSFER) 1,205,1 
205 ASSIGN 21 TO NSTART 

C STATEMENTS 7 TO 9 INITIALIZE NREV1 AND NREV2 FOR USE IN PART 7A. 

IF ( RB ( 2 ) ) 7,6,8 

6 IF ( VX( 2 ) ) 7,8,8 

7 ASSIGN 37 TO NREV1 
ASSIGN 35 TO NREV2 
GO TO 9 

8 ASSIGN 33 TO NREV1 
ASSIGN 37 TO NREV2 

9 DO 10 J=1,NEQ 
XDOTPM ( J , 1 ) = XDOTIJ) 

XINC(J) = 0. 

10 CONTINUE 

11 KSUB = 1 
ASSIGN 16 TO N 

C 

C PART 2. RUNGE-KUTTA SUBINTERVAL SCHEME. EQUATE PRODUCES THE NECCESSARY 

C DERIVATIVES XDOTIJ). 

12 DO 13 J=l,NfcQ 

XKIJ) * XDOTIJ) • DELT 

XINCIJ) = X INC ( J ) + AW(KSUB) *XK( J) 

13 XU) = XPRIM[J,2) + AK ( KSUB ) *XK ( J ) 

14 CALL EQUATE 

15 GO TO N, ( 16 , 17,18,20) 

C 

C PART 3. SUBINTERVALS 2, 3, AND 4, TO STATEMENT 19 FINISH A 

C RUNGE-KUTTA STEP AND INCREMENT XPRIM(J,2) IN DOUBLE PRECISION. 

16 KSUB = 2 
ASSIGN 17 TO N 
GO TO 12 

17 KSUB = 3 
ASSIGN 18 TO N 
GO TO 12 

18 DO 19 J=1,NEQ 

XINCIJ) = XINCIJ) ♦ AW! 4) *XDOT(J) * DELT 
180 CALL EXADDI XPRIMI J,2) , XPRIMB(J,2), XINCIJ)) 

XIJ) * XPR IM I J , 2 ) 

19 CONTINUE 
C 

C PART 4. BEGIN A NEW RUNGA-KUTTA STEP. THIS ALSO GIVES DERIVATIVES 

C FOR THE LOWER ORDER INTEGRATION CHECK. 

ASSIGN 20 TO N 
GO TO 14 

20 GO TO NSTART, (27,23,21) 

C 

C PART 5. STARTING PHASE PROGRAM. 

C PART 5 A . THIS SECTION COMPLETES THE FIRST STEP OF STARTING PHASE. 

21 ASSIGN 23 TO NSTART 
DO 22 J = 1 , NEQ 

OLD INC ( J ) =XINC I J ) 

XINCIJ) =0. 

XDOTPM I J , 2 ) = XDOTIJ) 

22 CONTINUE 
GO TO 11 


56 


C PART 5B. MAX ERROR TEST — STARTING ONLY — CHECK THE MAX ERROR AND 

C EITHER ENTER RUNNING MODE OR REPEAT START WITH SMALLER STEP. 

23 DO 2* J«2,NEQ 

24 XI NCI J ) ={X1NC( JJ+OLDINCI J ) ) *3.- ( XDOTPM ( J , 1 H-XDOTPM ( J , 2 > *4. 
l+XOOTI J) )*DtLT 

240 CALL ERROR/ 

25 IFIE2-ERLIMT) 26,26,56 

26 ASSIGN 27 TO NSTART 
ASSIGN LI TO I BEG IN 
A1 * A2 

GO TO 31 
C 

C PART 6. RUNNING PHASE PROGRAM. 

C PART 6A. CHfcCK THE INTEGRATION BY INTEGRATING OVER THE LAST 

C RUNGE KUTTA STEP BUT USE OOTS FOR LAST TWO INTERVALS, OLDDEL 

C AND CELT RESPECTIVELY. STATEMENT 28 IS THE LOWER INTEGRATION 

C MINUS RUNGE-KUTTA INCREMENTS. ERRORZ COMPUTES THE MAXIMUM RELATIVE 

C ERROR AND STATEMENT 29 TESTS THIS ERROR AGAINST THE LIMIT VALUE. 

27 RATIO * DELT/OLODEL 
HFACT-DELT/I 1.+RAT10) 

ACGEF1«— RAT 1 0«RAT IQ*HFACT 
AC0EF2*RAT 10* ( DEL T+3. *OLDDEL ) 

AC0EF3-DELT+0ELT+HFACT 

DO 28 J=2, NEQ 

28 XINC(J) = ACOEF1*XOOTPMIJ,1) +AC0EF2*XDDTPMt J,2)-6.*XINCIJ) 

1+AC0EF3»X0QT I J ) 

280 CALL ERRORZ 

29 IF IE2-ERLIMT) 30,30,57 
C 

C PART 7A. LAST POINT OKAY. COUNT THE REVOLUTIONS PAST THE X-AXIS. 

C A STEP GREATER THAN 1/2 REV. MAY FAIL TO ADO IN. 

30 H2 = OELT 

31 QMAX = MAX IF ( Q, QMAX ) 

IF 1 RBI 21 ) 32,34,34 

32 GO TO NREV1, (37,33) 

33 ASSIGN 37 TO NREV1 
ASSIGN 35 TU NREV2 
GO TO 37 

34 GO TO NREV2 , 137,35) 

35 ASSIGN 33 TO NREV1 
ASSIGN 37 TO NREV2 

36 REVS = REVS + 1. 

37 IF ( XABSF( IM0DE)-1) 42,38,42 
C 

C PART 7B. IN ORBIT ELEMENTS. ADJUST ARGUMENT OF PERICENTER AND MEAN ANOMALY 

C TO + OR - PI TO MAINTAIN ACCURACY IN SIN-COS ROUTINES. 

38 IF (EMONE) 39,42,42 

39 DO 41 J=4,7,3 

ADJ2=INTF( XPRIM( J , 2 ) /6. 283 18532+S I GNF ( .5,XPRIM( J,2» ) ) 

IF IADJ2) 40,41,40 

40 ADJ3 = -ADJ2*6. 28125 

400 CALL EXAOD( XPRIMI J, 2) , XPRIMB( J,2) , ADJ3) 

AD J3=-ADJ2*. 00 19353072 

401 CALL EXADDUPRIMl J , 2 ) , XPR I MB ( J , 2 ) , AD J3 ) 

41 CONTINUE 
C 

C PART 7C. ADVANCE THE REMAINING PARAMETERS, FIND NEW STEP SIZE, 

C AND TEST FOR AN ORIGIN TRANSLATION. 

42 DO 43 K= 1 , 3 
XWHOLE(K)=VX(K) 

43 XWHOLE ( K+3 ) = RB(K) 

DO 44 J- 1 1 NEQ 

XDOTPM I J , 1 ) = XDOT PM ( J , 2 ) 

XDOTPM I J , 2 ) = XDOTIJ) 

XPRIMI J , 1 ) = XPRIMI J,2) 

XPR I MB ( J , 1 ) = XPRIMBl J,2) 

XINCIJ) = 0. 

44 CONTINUE 
OLDDEL * DELT 

45 CALL STEP 

IF (DONE) 67,450,67 

450 IF (NSTART) 451,1,451 

451 IF ( MBQDYS ) 46,47,46 

46 CALL TESTTR 

IF ( TRSFER ) 1,47,1 

47 IF ( XABSFl IMODEI-3) 11,48,11 
C 

C PART 7D. IF IN TEMPORARY RECTANGULAR COORDINATES, TEST FOR RETURN 

C TO ORBIT ELEMENTS. FIRST, E IS FOUND. IF TIME HAS NOT ADVANCED 

C SUFFICIENTLY, INTEGRATION CONTINUES IN RECTANGULAR VARIABLES ( STATE. 48). 

C STATEMENT 49 DETERMINES IF KEPLERS EQUATION CAUSED I MODE = 3. IF NOT, 

C AN E CLDSE TO 1 CHECK IS MADE IN STATEMENT 50. IF IT DID, RECTANGULAR 

C VARIABLES WILL BE USED IF THE LIMIT IS TOO SMALL (STATEMENT 52), OR 

C IF E IS 5 OR GREATER (STATEMENT 53) OR IF THE PATH LIES CLOSE TO AN 

C ASYMPTOTE (STATEMENT 55). 

48 CALL CON VT l (VX.AMC) 

EXMODE=SQRTF ( 1 . +AMSQRD/GK2M* ( VSQRD/GK2M-2. /R ) ) 

EMONE=EXMODE— l« 

IF ( (XPRIMI 1 l-TTEST) *DELT) 11,11,49 

49 IF (ASYMPT) 51,50,51 

50 IF ( ETOL-ABSFl EMONE) ) 55,11,11 

51 I F ( EMONE ) 55,55,52 

52 IF(CONSTU— 1. E— 7 ) 11,53,53 

53 IF (EXMQDE-5.) 54,11,11 

54 CALL C0NVT2 

IF ( ABSF ( TRU )— 2.2/ SQRTF ( EXMODE ) ) 55,55,11 

55 ASYMPT = 0.0 
IM0DE*~2 

555 CALL TESTTR 
GO TO l 
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PART 8. COMtS HERE WHEN ERROR TEST FAILED — BOTH STARTING AND RUN. 

RETRIEVE OLU POINT AND RECOMPUTE WITH SMALLER INTERVAL. 

IF TWO CONSECUTIVE TRYS FAIL (STATEMENT 59) THE STARTING SEQUENCE OCCURS. 

56 ASSIGN 1 TO IBEGIN 

57 DO 58 J=1,NEQ 

XPR IH ( J » 2 ) = XPRIM(J,1) 

XPRIMB(J,2) = XPR I MB ( J » l ) 

XDOTIJ)=XDOTPM( J,2) 

X INC ( J ) = 0. 

58 CONTINUE 
STEPN0=STEPN0+1. 

H2 = DELT 

DELT=S IGNF ( EXPF ( ( ERL0G-A2 ) /5. ) » DELT ) 

A2 =A I 

59 IF (FAIL-STEPGO) 60,61,60 

60 FAIL = STEPGO 

GO TO IBEGIN, (11,1) 

61 ASSIGN 1 TO IBEGIN 

IF ( STEPNO + STEPGO - STEPMX) 62,62,45 

62 GO TO IBEGIN, (11,1) 

PART 10. PRINT OUT THE ERROR INFO. IF EREF HAS A - SIGN. THEN RETURN. 

67 IF (EREF) 60,72,7 2 

68 WRITE OUTPUT TAPE 6,70 
REWIND 4 

DO 69 1=1, INDERR 
READ TAPE 4, BEX 

69 WRITE OUTPUT TAPE 6, 71, BEX 
REWIND 4 

INDERR = 0 

70 FORMAT ( 7H1 STEP,6X,4HT IME»6X»4HDELT,7X» 2HA2 » 8X , 2HE2 , 7X, 4HMASS » 6X» 

14HE, VX,4X, 8H0M£GA»VY,2X»8HN0DES»VZ»3X,6HINCL»X,5X»4HMA, Y »6X»3HP, Z , 

24X , 1HK// ) 

71 FORMAT ( F5. , 1H+F3. , 1P11G10.2, 12 > 

72 RETURN 
ENO 

SUBROUTINE ORDER 

THIS ROUTINE TAKES THE BODY LIST READ FROM CARDS AND SORTS THEM IN 
ORDER SO THAT THE DISTANCE FROM THE REFERENCE TO EACH BODY IS 
DEPENDENT UPON ALREADY COMPUTED DISTANCES ONLY. 

ELLIPSE DATA ARE READ INTO A BLOCK OF 120 STORES RESERVED FOR 
TEN ELLIPSES. ONE ELLIPSE IS READ INTO A 12 STORE BLOCK. 

THE SINES AND COSINES OF THE 3 ANGLES ARE COMPUTED AND STORED 
IN THE TDATA ARRAY ALONG WITH THE REST OF THE ELLIPSE DATA. 

A BLOCK IS ARRANGED AS FOLLOWS' 

(1) = NAME OF BODY IN BCD, ONLY 6 CHARACTERS. 

(2) = NAME UF REFERENCE BODY IN BCD, SAME RESTRICTION. 

(3) = MASS OF THE BODY IN SUN MASS UNITS. 

(4) = RADUIS INSIDE OF WHICH COORDINATES WILL BE TRANSLATED TO THIS BODY. 

(5) = SEMILATUS RECTUM IN ASTRONOMICAL UNITS. 

(6) = ECCENTRICITY OF THE ORBIT. 

(8) = LONGITUDE OF ASCENDING NODE. 

(7) = ARGUMtNT OF PERIHELION. 

(9) = INCLINATION OF THE ORBIT. 

(10) = PERIGEE PASSAGE JULIAN DAY. 

(11>= PERIGEE PASSAGE FRACTION OF DAY. 

(12)= PERIOD OF THE ELLIPSE IN MEAN SOLAR DAYS. 

AMASS = MASS OF EACH BODY, SUN MASSES. ORDER OF PNAHE. 

BMASS = SELECTED FROM AMASS. CORRESPONDS TO BNAME LIST. 

BNAME = THE ORDERED LIST OF BCD BODY NAMES. CAN BE USED IN OUTPUT . COMMON. 
BOOYCD = THE ORIGINAL BCD NAMES READ FROM CARDS. 

BODY L = THE LIST OF BCD BODY NAMES WITH THE REFERENCE BODY AT TOP. 

INITIALLY EQUAL TO BODY CARD LIST (BOOYCD). 

I BODY = ARRAY OF SUBSCRIPTS. WHEN A DISTANCE IS FOUND FROM EPHEMERI S , IT 
MAY BE ADDED (OR SUBTRACTED) FROM THE BODY POSITION GIVEN BY 
XP(IBODY) TO OBTAIN THE POSITION OF THE PRESENT BODY. COMMON. 
KZERO = COUNT OF ZERO REFERENCES. THERE MUST BE ONE AND ONLY ONE ZERO. 

FROM LOCATION IN BNAME LIST. NOT IN COMMON. 

MANE = ARRAY OF SUBSCRIPTS. INVERSE OF NAME. GIVES NEW LOCATION OF 
BNAME LIST IN TERMS OF BODYL. NOT IN COMMON. 

NBQDYS = COUNTED INTERNALY. TOTAL NUMBER OF BODYS. 

MBOOYS = COMPUTED INTERNALY. TOTAL NUMBER OF EPHEMER IDES (NBODYS-l). 

NAME = ARRAY OF SUBSCRIPTS. GIVES OLD LOCATION OF NAMES IN BODYL 
NEFMRS = ARRAY OF SUBSCRIPTS. GIVES LOCATION OF BODY IN PNAME LIST 
IN TERMS OF THE EFMRS LIST. STORED IN COMMON. 

NREFER = ARRAY OF SUBSCRIPTS. LOCATES THE REFERENCE BODY IN BODYL. 

ORDER OF THE ARRAY CORRESPONDS TO BODYL. NOT IN COMMON. 

NNREFR = ARRAY OF SUBSCRIPTS. LIKE NREFER BUT REFERS AND CORRESPONDS TO 
BNAME LIST. NOT IN COMMON. 

PNAME = A PERMANENT LIST OF BCD BODY NAMES. 1 WORD EACH 16 CHARACTERS 

MAX). USED TO IDENTIFY MASS, REFERENCE NAMES, ETC. THE LIST IS 
A MAXIMUM OF 30 NAMES. PRECISION TAPE NAMES ARE FROM l TO 20, 
ELLIPTIC NAMES ARE FROM 21 TO 30. 

REFER = A PERMANENT LIST OF BCD BODYS THAT ARE THE REFERENCES OF 

DISTANCES GIVEN IN EPHERMER IDES (TAPES OR ELLIPSE). CORRESPONDS 
TO PNAME LIST. 

COMMON C 


DIMENSION 

A ( 600 ) , 

B ( 700 ) , C( 4000 ) , 



1 

AMASS 

(30) , 

BMASS 

(8) , 

BNAME 

(8) , 

2 

BODYL 

(8) , 

EFMRS 

(7), 

IBOOY 

(8) , 

3 

MANE 

(8) , 

NAME 

(8) , 

NEFMRS 

(8) , 

3 

NEFMRT 

(8) , 

NNREFR 

(8) , 

BOOYCD 

(8) , 

4 

NREFER 

(B), 

PNAME 

(30) , 

RBCRIT 

(7), 

5 

RCRIT 

(30) , 

REFER 

(30) , 

TDATA 

( 18,7) 

6 

TDEL 

(7) , 

TIM 

(7), 

EL I PS 

( 120) , 

7 

NDUD 

(9), 

XPR I M 

(200) 
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EQUIVALENCE 


1 ( A 

,C 

I 111), (AMASS , A 

( 347 ) ) , ( A TMN ,A 

( 21)), 

2 ( AU 

, A 

( 29)), IB ,C 

( LLLL 1 ) t (BMASS ,B 

( 1371), 

3( BNAME 

,B 

( 122) ) , ( BODYCD,A 

l 143 ) ) , ( BOOYL ,B 

( 153)), 

4 ( E FMRS 

,8 

( 130) >, (BLIPS , A 

( 167)), (FILE ,B 

( 22)), 

5 ( GK2M 

, B 

( 36 ) ) , ( GKM , B 

< 37) ) , ( I BODY , B 

( 177)), 

6 (MBODYS 

, B 

( 42) ), (NBOOYS, B 

( 4 l ) ) , 1 NEFMRS , B 

1 185)), 

7 ( PNAME 

, A 

( 287) ) , (RATMOS.B 

( 23) ) , ( RATM , A 

( 22)), 

8IRBCRIT 

, B 

( 1 45 ) ) , ( RCR I T , A 

( 377)), (REFER ,A 

( 317)), 

91RE , A 

EQUIVALENCE 

( 25) ) , (RfcSQRD.B 

( 7) ), (REVOLV.B 

( 21) ) 

KROTATE 

, A 

( 39) ) , ( SPD , A 

( 44) ) , ( SQR0K1 , A 

( 47)), 

2 ( SQRDK 

,0 

( 35) ) , ( TDATA ,B 

( 265) ), ( TDEL ,B 

( 170)), 

3ITFILE 

» A 

{ 6)), (TIM , B 

( 163) ), (TRSFER, B 

( 8) ) , 

4UPRIM 

,C 

( 711) ) , (OUTPOT.B 

( 399)) 



EQUIVALENCE t MANE ( i J . NOUDI 2 U 

THIS SECTION SEES WHAT ELLIPSE DATA WAS READ FROM CARDS AND PUTS THE 
NAMES IN PLACE SO THAT DATA WILL BE USED IF NEEDED. ELLIPSE DATA HAS 
PRIORITY OVER TAPE DATA BECAUSE LAST DATA IN LIST IS THAT ACTUALLY USED. 
FUNCTION COMP ARF ( A * B ) IS EQUIVALENT TO <A-B) BUT WILL NOT OVERFLOW. 

COMP ARF ( A , B ) = ( A+B ) * I - ( A* B ) ) 

DO 2 K= 1 * 120 1 12 
I F ( EL I PS ( K ) ) 1,2,1 
KOUNT * ( K— 1 ) / 12+2 l 


PNAME I KOUNT ) 
REFER ( KOUNT ) 
AMASS(KOUNT) 
RCR I T ( KOUNT ) 
CONTINUE 

PART 0. 


EL IPSIK) 
ELIPS(K+1) 
ELIPSIK+2) 
ELlPSIK+3) 


THROW AWAY BLANKS AND DUPLICATES IN BNAME LIST. 

ALSO COUNT THE BODIES. 

IF (TRSFER) A, 3,4 

3 BNAME(l) - BOOYCD(l) 

4 DO 5 K= 1 i B 

5 BNAME(K+1)= BODYCO(K) 

L = l 

BOOYL(O) = 0. 

DO 8 1=1,9 
BODYLt I) = O. 

DO 6 K= 1 , L 

IF (COMPARF { BNAME ( I ) » BODYL ( K— 1 ) } ) 6,7,6 

6 CONTINUE 

BODYL ( L ) = BNAME ( I > 

L = L+l 

7 BNAME ( I) = 0. 

8 CONTINUE 
NBODYS = L— 1 
MBODYS = NBOOYS-1 

PART 1. FIND THE REFERENCE BODY FOR EACH BODY IN THE LIST OF BODYS 
READ FROM CARDS. CLEAR NREFER ANO BNAME. 

DO 13 KL=1, NBODYS 
NREFER(KL) = 0 
NEFMRTIKLI =0 
BNAME CKL) = 0. 

DO 12 KP= 1,30 

IF I COMPARE! BODYL! KL ) » PNAME ( KP ) ) ) 12,9,12 

9 NEFMRT(KL) = KP 

DO 11 KK = 1,8 

IF ( COMPARE! REFER! KP) »BQDYL(KR))) 11,10,11 
10 NREFER(KL) = KR 
1.1 CONTINUE 

12 CONTINUE 

13 CONTINUE 

PART 2 . COUNTS 0 REFERENCES AND SAVES TEMPORARY SET OF INDEXS. 

14 IF (NBOOYS) 24,24,15 

15 KZEROS = 0 
MISPEL = 0 

DO 20 K = 1, NBODYS 
NNREFR(K) = NREFER ( K ) 

16 IF ( NEFMRT( K ) ) 18,17,18 

17 MISPEL = MISPEL + 1 

18 I F ( NREFER IK)) 20,19,20 

19 KZEROS = KZEROS ♦ 1 

20 CONTINUE 

21 IF I KZEROS- 1) 24,22,24 

22 IF (MISPEL! 24,23,24 

23 IF I NBODYS— 8 ) 28,28,24 

PART 3 . REPORTS ERRORS IN BODY LIST. 

24 WRITE OUTPUT TAPE 6,25 , NBODYS ,M I SP6L , KZ EROS ,( BODYL ( K ), K= 1 , NBODYS ) 
WRITE OUTPUT TAPE 6,26 , I NREFER < K ) , K= 1 , NBODYS ) 

WRITE OUTPUT TAPE 6,27 , { K , PNAME I K ) , REFER I K > , K=1 , 30 ) 

25 FORMAT (26H0G00FY BODY LIST (NBODYS =I2,13H, MISSPELL =12, 

1 11H, KZEROS =12, lHI/ii HO BODY LIST =8(3X,A6M 

26 FORMAT C11H NREFER =16,719) 

27 FORMAT I/5I3H K3X, 4HB0DY4X, 5HREF6R5X, 1/51 I3,2X,A6,2X,A6,5X) ) 

CALL EXIT 

PART 4. TRACES OUT ..REFERENCE TO BODY.. RELATIONSHIPS 

28 KK = 2 
KN * 1 
NAME ( 1 ) = 1 

29 IF ( NREFER ( KN ) ) 24,31,30 

30 NAME IKK) = NNREFR ( KN ) 

NNREFR(KN) = 0 

KN = NAME ( KK ) 

KK * KK ♦ 1 
GO TO 29 
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C PART 5. TRACES OUT ..BODY TO REFERENCE.. RELATIONSHIP 

31 DO 34 KN = 1 • NBODYS 
DO 34 K = 1 * NBODYS 

32 IF (NNREFR(K) - NAMEIKN)) 34,33*34 

33 NAME IKK) = K 
KK = KK + 1 

34 CONTINUE 
C 

C PART 6. INVERTS NAME TO MANE, STORES BNAME, BMASS * RBCRIT, AND A 

C TEMPORARY NEFMRS. 

DO 35 K = 1 * NBOOYS 
N = NAME IK) 

MANE(N) = K 
NEF * NEFMRT(N) 

BNAME ( K ) = PNAME(NEF) 

BMASS ( K ) = AMASS ( NEF ) 

RBCRIT(K) = RCR I T ( NEF ) 

NEFMRS ( K ) » NEF 

35 CONTINUE 
C 

C PART 7. FINDS NNREFR REFERENCE FOR BNAME LIST , ALSO TEMP. IBODY 

DO 36 K = 1, NBODYS 
N * NAME! K ) 

NRF = NREFER(N) 

NNREFR ( K ) = MANE ( NRF ) 

36 I BODY! K ) = MANE ( NRF J 
C 

C PART 8 . FINDS IBODY FOR BACKWARD REFERENCE. 

DO 39 K=l, 8 

37 IF ( NNREFR! K ) ) 24,40,38 

38 N = NNREFR ( K ) 

I BODY ( N ) * -K 

39 CONTINUE 

C IBODY LIST IS COMPLETE. 

C 

C PART 9 . WRITES OUT EPHEMER I S LIST TO BE USED IN STORING DATA AND 

C MAKES FINAL NEFMRS LIST. 

40 KK = 1 

DO 43 K=l, NBODYS 

41 I F ( NNREFR ( K ) ) 42,43,42 

42 EFMRS(KK) = BNAME t K ) 

NEFMRS ( KK ) = NEFMRS ( K ) 

KK = KK + 1 

43 CONTINUE 

NEFMRS ( NBODYS ) = 0 
C 

C PART 10. SAVES ELLIPSE DATA 

FILE = 0. 

IF IMBODYSi 430,480,430 
430 DO 48 K=1 , MbGDYS 

44 I F 1 NEFMRS ( K )— 20 ) 47,47,45 

45 DO 46 J=5, 12 

L= ( NEFMRS ( K ) - 21) • 12 +J 

46 TDATAIJ— 4,K) = ELIPS(L) 

DO 50 J=7, 9 

L = (NEFMRS I K )— 2 1 ) *12+ J 
TDATA ( J + 2, K ) = S INF ( EL I PS( L ) ) 

50 TDATA( J+5 » K ) = COSFI EL IPS( L ) ) 

GO TO 48 

C 

C PART 10A. LOADS A FALSE (VERY EARLY) TAPE TIME TO FORCE TAPE 

C READING BY THE EPHMRS ROUTINE. FILE = 0 UNLESS TAPE IS USED. 

47 TDEL(K) = 0 
TIM(K) = 2400000.5 
FILE = 10. 

48 CONTINUE 
C 

C PART 11. COMPUTE GRAVITATIONAL CONSTANTS. 1.9866 E+30 = KILOGRAMS/SUN MASS 

C IF ORIGIN BODY HAS AN ATMOSPHERE, SET ROTATION RATE AND ATMOSPHERE RADIUS. 

C POSITION THE EPHEMER IDES TAPE AT THE BEGINNING OF THE CORRECT EPHEMER I S 

C BY MATCHING THE EPHEMER I S NUMBER READ FROM TAPE (FILE) WITH THE DESIRED 

C EPHEMER I S NUMBER 1TFILE). 

C 

480 RESQRD = RE**2 

SQRDK = SQRUK1*AU**3/SPD**2 

GK2M = SQRDK* ( BMASS ( 1 ) + XPR IM I 2 ) / l .9866 E30) 

GKM = SQRTF ( GK2M ) 

RE VOL V = 0. 

RATMOS = 0. 

IF ( ATMN-BNAME (1)1 51,49,51 

49 REVOLV = ROTATE 
RATMOS = RATM 

51 IF (FILE) 56,56,52 

52 CALL BSFILEI3) 

53 READ TAPE 3, FILE 

IF (FILE-TFILE) 54,56,55 

54 CALL SKFILEI3) 

GO TO 53 

55 BACKSPACE 3 
BACKSPACE 3 
GO TO 52 

C 

C PART 12. WRITES THE BNAME LIST ON TAPE 6 . 

56 IF (OUTPOT) 58,59,58 

59 WRITE OUTPUT TAPE 6,57, BNAME ( 1 ) , ( BNAME 1 K ) » K«2 » NBODYS ) 

57 FORMAT ( 19H0REFERENCE BODY IS A6,5X,23H PERTURBING BODIES ARE 

1 7 ( 2X , A6 ) ) 

58 RETURN 
END 
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subroutine oblate 

C THIS SUBROUTINE COMPUTES THE OBLATENESS ACCELERATIONS (OBLAT) DUE TO AN 
C AXIALLY SYMMETRIC EARTH. THE 2ND, 3RD, AND 4TH SPHERICAL HARMONIC COEFF. 
C ARE OBLATJ » OBLATH, AND OBLATD RESPECTIVELY. 

C 

COMMON C 
C 

DIMENSION A ( 600) , BI700), C(4000>, 

1 RBI3), OBLAT ( 31 

C 

EQUIVALENCE 

1 I A ,C ( II)), (B ,C C I 111 ) ) , { GK2H ,8 ( 361), 

2 I OBLAT J, A ( 26) ) 1 1 OBLATD, A ( 27 ) ) , ( OBLATH, A l 28)), 

3 ( OBLAT ,B I 75)), (R ,B ( 1021), (RB ,B I 193)), 

4 ( RE , A ( 25) ) , (RSQRD ,B ( 45 ) ) , ( R6SQRD, B ( 7)) 

C 

AA = RB(3)/R 
AB = AA* AA 

IF ( ABSF ( AA)— 1 . E— 6) 1.1,2 

1 AA=0. 

AB=0. 

2 AC * RE SQRO/ RSQRD 

AD - GK2M/RSQRD/R*AC 
AE * OBLAT J*AD 
AF * OBLATH*AD*RE/R 
AG - OBLATD*AD*AC 

AH * AE* < 5. AB-l. )+AF*( 7. AB-3. ) *AA+AG*( 6. AB-9.AB* *2-0. 4285714286) 
OBLAT(l) * AH*RB 1 1 ) 

OBLAT ( 2 ) - AH*RB12) 

OBLAT ( 3 ) = ( AH— 2.AE+AG* ( 4.AB— 1. 714285714) )*RB( 3)-AF* ( 3. AB-0.6) *R 

3 RETURN 
END 


SUBROUTINE OUTPUT 

C ENTS AND RECTANGULAR COORDINATES ARE OUTPUTTED. IF THE OBJECT IS NOT WITH 

C THIS IS THE ROUTINE WHICH FORMS THE BASIC DATA OUTPUT. BOTH ORBIT ELEM- 

C IN AN ATMOSPHERE (PRESS=0.), ONE LINE OF DATA IS DELETED. LIKEWISE, 

C ONLY THOSE PERTURBING BOOIES PRESENT HAVE THEIR DISTANCES OUTPUTTED. 

C 

COMMON C 

C 


DIMENSION 

A ( 600 ) 

, B ( 700 ) , C14000), 



1 R 

(8) , 

ORBELS (6), 

VATM (3), 


2 BNAME 

(8) , 

RB ( 3,8), 

DIRCOSl 3,8) 


3 XPR I M 

(200) , 

RAMC (5) 



EQUIVALENCE 




11A ,C 

l 

ID), (ALPHA , A 

( 49) ) , (ALT , A 

( 4) ) 

2 ( AMC ,B 

( 

87)). CAM ,B 

( 901), (AREA ,B 

( 6) ) 

3 ( BNAME , B 

( 

122)), IB ,C 

1 1111) ) , (CD , A 

( 165)) 

4 (CL , A 

( 

164 ) ) , ( COSALF, B 

( 48) ) , (COSTRU, B 

( 53) ) 

5 ( DTOFF J , A 

1 

23)1, (H2 , B 

( 15) > , I IMODE , A 

( 1) ) 

6 ( MBODYS , B 

l 

42) ) , ( NBODYS, B 

< 41) ) , (ORBELS, B 

( 116)) 

7 ( PRESS ,B 

( 

33)), IP , B 

( 84)), (PSI , B 

( 30) ) 

8 ( PSI R ,B 

I 

398) ) , (PUSH , A 

( 166)), (Q , B 

( 59) ) 

9{ RAMC ,B 

( 

393 ) ) * ( RB ,B 

( 193)), (REVS , A 

t 48) ) 

EQUIVALENCE 




1 ( R , B 

( 

102) ) , (SINALF.B 

( 46) ) , (SINTRU.B 

( 52) ) 

2 ( STEPGO , A 

t 

41)), (STEPNO, A 

( 42 ) ) , ( TABLT ,B 

( 20) ) 

3ITRU ,8 

( 

40)), (VATM ,B 

( 97) ) , (VQ ,B 

( 100)) 

4 ( V ,B 

( 

95 ) ) t ( VX , B 

( 92)), (VY , B 

l 93) ) 

5 ( VZ , B 

( 

94) > , (XPRIM ,C 

( 711)), (OUTPOT, B 

( 399)) 


c 

DAYJ= 1 DTOFF J— 2.4E6 ) +TABLT 
ALPHA 1 = ALPHA*57. 29577951 

REV = REVS+ARCTANl— RBI2) ,-RB( 1) )/6. 28318532 + .5 

16 CALL CONVTK VX, AMC ) 

IMODE=IMODE 

GO TO (2,1,1), IMOOE 

1 CQ0E=6HRECTAN 

18 CALL CONVT 2 
GO TO 4 

2 DO 3 K=l,6 

3 ORBELS(K) = XPRIMIK+2) 

C0DE=5H0RB I T 

TRU=ARCTAN( S INTRU, COSTRU ) 

4 PSI = ATANF1 (RB( 1 ) * VX+R6 ( 2 ) * VY+RB ( 3 ) • VZ ) / AM ) 57. 2957795 
IF (OUTPOT) 19,6,19 

6 WRITE OUTPUT TAPE 6, 11 , STEPGO, STEPNO, ORBELS ( 1 ) , ORBELS ( 2 ) , V, R 1 1 ) , B 
I NAME ( 1 ) , CODE , I MODE ,XPRIM(1),0RBELS(6) ,TRU,VX,RB( 1 ) , XPR I M ( 2 ) , DA Y J , 0 
2RBELS l 5) .ORBELS! 3) , VY , RBI 2 ) ,REV , ALPHA 1 , PSI ,QRBELSI 4) , V2 »R8 13) , H2 
C 

C IF WITHIN AN ATMOSPHERE COMPUTE DRAG, LIFT, G, ETC., AND PRINT EXTRA LINE. 

19 IF (PRESS) 5,7,5 

5 X I FT = Q*AREA*CL 
ORAG = Q*AREA*CD 

G “ (PUSH-0RAG*C0SALF+XIFT*SINALF)/XPRIM(2)/9. 80665 

17 CALL CONVT 1 ( VATM, RAMC ) 

PS 1R * ATANF( (RB(U*VATM(U+RB(2)*VATM(2)+RB(3)*VATMt3) )/RAMC(4) )* 

1 57.2957795 
IF (OUTPOT) 7,14,7 

14 WRITE OUTPUT TAPE 6, 12 , ALT , PS IR, DRAG, VQ , G, PUSH 
C 
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C IF PERTURBATING BODIES ARE PRESENT, FIND THEIR DISTANCES AND PRINT THEM. 

7 IF { MBOD YS) 8,10,8 

8 DO 9 J=2 , NBODYS 
DO 9 K= 1 , 3 

9 D I RCOS ( K , J ) = -RBIK, J)/R( J) 

IF (OUTPOT) 10,15,10 

15 WRITE OUTPUT TAPE 6,13, 

llBNAMEI J),R( J),OIRCOS( 1 , J ) , D I RCOS ( 2, J ) , D I RCOS ( 3 , J ) , J=2 , NBODYS ) 

10 RETURN 

11 FORMAT ( 6H0STEP=F6« , 2H +F4. , 3X, 13HECCENTR IC I T Y= IPG15 . 8 , 7H QMEGA=G 15 

1.8.4H V=G15. 8 , 3H R=G15.8,7H REFER=A6, IX, A6, I2/6H TIME=1PG14.7, 14 

2H SEMILATUS R.=G15.8,7H TRU A=G15.8,4H VX=G15.8,3H X=G15.8,7H RMAS 
3S=G15. 8/9H JDAY= 240PF10.4, 15H MEAN ANOMALY 3 IPG 15. 8 , 7H NODE=G15. 

48, 4H VY=G15. 8, 3H Y=G15.8,7H RE VS. -G15. 8/ 6H ALF A=G14 . 7 , 14H PATH A 
5NGL6=G15.8,7H INCL=G15.8,4H VZ=G15.8,3H Z=G15.8,7H DELT=G15.8) 

12 FORMAT ( 6H ALT. = 1PG14. 7, 14H R PATH ANGLE=G1 5. 8, 7H DRAG=G15. 8 , 4H VR 

1=G15.8,3H G=G15.8,7H PUSH=G15.8) 

13 FORMAT (2(1X,A6,3H R= l PG1 4. 7 , 0P3F 10 . 6, 1 IX ) ) 

END 

FUNCTION QUAD (X,IC) 

C THIS ROUTINE COMPUTES ANY VARIABLE, QUAD, AS A QUADRATIC FUNCTION OF X. 

C QUAD = A + BX + CXX. THERE MAY BE SEVERAL SETS OF COEFF IENTS, EACH SET 

C BELONGING TU A PARTICULAR REGION OF X. THE COEFN ARRAY IS ARRANGED AS — 

C XI, A1,B1,C1,X2,A2,B2,C2,X3,A3,B3,C3,X4, 

C WHERE A 1 , 81 , Cl ARE THE COEFF IENTS TO BE USED FOR X BETWEEN XI AND X2,ETC. 

C ANO XI IS LESS THAN X2, X2 IS LESS THAN X3, X3 IS LESS THAN X4, ETC. 

C IC IDENTIFIES WHICH DEPENDENT VARIABLE, QUAD, IS BEING SOUGHT. 

C ICC(IC) DEFINE THE STARTING LOCATIONS IN THE COEFN ARRAY FOR VARIABLES X. 

COMMON C 
C 

DIMBNSION A ( 600 > , 81 700) , Cl 4000), 
l COEFN ( 19U ) , ICC ( 5 ) 

C 

EQUIVALENCE 

1 ( A ,C ( 1 1 ) ) , ( B ,C (1111)), (COEFN ,A ( 407)), 

2 ( ICC , A ( 153)) 

C 

I=ICC( IC) 

1 IF (X-COEFN(I)) 2,3,3 

2 I - 1-4 
GO TO 1 

3 I F ( X— COEFN ( i +4 ) ) 5,5,4 
41=1+4 

GO TO 3 

5 QUAD = COEFN ( I+1)+X*(C0EFN( 1+2 )+ X* COEFN ( 1+3) ) 

ICC C IC) = I 

RETURN 

END 


SUBROUTINE STAGE 

C THIS ROUTINE IS CALLEO TO PREPARE DATA FOR USE IN NBODY. STAGE DATA IS 

C TAKEN FROM PERM INENT STORES ANO LOADED INTO WORKING STORES. STAGE DATA 

C MAY BE SET ASIDE FOR LATER USE (IF ON NSAVE-NSTAGE ) . WHEN IMODE IS 4, 

C CONVERSION FROM EARTH-SPHERICAL TO RECTANGULAR OR ORBIT ELEMENTS TAKES 

C PLACE IN TUBES. 

C 

COMMON C 
C 

DIMENSION A { 60 0 ) , BI700), C(4000), 

1 XPR I M ( 200 ) ,XPRIMB(200) ,TB( 10 ) , FLOW 1 ( 10 ) , AEX I T II 1 0 ) , S IMP 1 ( 10 ) , 

2AREAK 10) ,DELT1 ( 10) , I DENT ( 10 ) , TA BLE ( 200 ) , RMASS 1 ( 10) , 0(600 I 
C 


EQUIVALENCE 


1 (A 

,c 

( ID), ( AEX I T ,B 

( 3) ) , (AEX IT1, A 

( 103)) 

2 ( AREA1 

i A 

< 113)), (AREA , B 

( 6) ) , ( B ,C 

( 1111) ) 

3 ( DELT1 

A 

( 133 ) ) , ( DELT , B 

( 1 ) ) , ( D ,C 

(2111) ) 

4 ( DEL 

• A 

( 43) ) , ( DELMAX , A 

( 19)), (DONE ,B 

( 39) ) 

5 ( EREF 

• A 

( 13 ) ) , ( ERLOG , B 

( 17) ) , (EXITA , B 

< 392)) 

6 ( FLOW 

B 

( 5 ) ) , ( FL0W1 , A 

( 83) ), ( IDENT , A 

{ 123)) 

71 IMODE 

A 

( 1 ) ) , ( LSTAGE, A 

1 38) ) , ( MODOUT , A 

1 20) ) 

8 ( NCASE 

r C 

( l) ) , INCASES, A 

( 600 ) ) , ( NSA VE ,C 

( 4) ) 

9 ( NSTAGE , A 
EQUIVALENCE 

I 3 ) ) , ( PUSHO ,B 

( 391) ) , (RMASS1,A 

( 73) ) 

11SIMP1 

A 

( 93)), (SIMP ,B 

( 2) ) , ( TB , A 

( 63) ) 

2 ( TABLE 

C 

(1911) ), (TKICK , A 

( 15)), UMAX ,B 

( 4) ) 

3( TTOL 

A 

( 45 ) ) , l XPR I MB ,C 

( 9 1 1 ) ) , ( XPR I M ,C 

( 711)) 

4 (RETURN 

6 

( 400) ), (OUTPOT, B 

( 3991) 



C 
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load stage data into forking 


i 


C PART 0. SAVE INITIAL DATA IF DESIREO. 

C STORAGE. ALLOW ADDITIONAL STAGE INPUT. 

IF (DEL) 100,99,100 
99 DEL * DELMAX-TKICK 

100 IF (NSAVE-NSTAGE) 103,101,103 

101 NCASES = NCASE 

00 102 J=l,liOO 

102 D( J ) * A ( J ) 

IF (OUTPOT) 103,97,103 

97 WRITE OUTPUT TAPE 6 , 98 .NSTAGE , NCASE 

98 FORMAT ( 29H ^AVED INITIAL DATA FOR STAGEI2 , 8H OF CASEI4.1H.) 

103 NSTAGE - NS I AGE 

TMA X = XPRIMt 1 )+TB(NSTAG£) 

XPRIMBI2) = 0. 

IF ( RMASS l (NSTAGE 1 I 117,117,118 

117 XPR I M( 2 ) = XPRIM121+RMASSUNSTAGE) 

GO TO 119 

118 XPR I M t 2 ) * RMASS1 ( NST AGE ) 

119 FLOW = FLOWUNSTAGEI 
SIMP * SIMP1 I NSTAGE ) 

AEXIT * AE X 1 T 1 1 NSTAGE ) 

AREA * AREAHNSTAGE) 

OELT * DEL T 1 1 NSTAGE } 

ID * I OENT ( NSTAGE I 

3 CALL INPUT { ID, C, TABLE) 

ERLOG * L0GF(A8SF( EREF) ) 

TTOL * 5E-8*A8SF l TMAX) ♦ 1 E— 8 
PUSHO » SIMP*FLQW*9. 80665 
EXITA = AEXIT* 100. 

MODGUT * MQUOUT 

IF ( DELT I 105,104,105 

104 DELT = TBlNiTAGEJ/lOO. 

OELT l ( NSTAGfc I = DELT 

105 GO TO 1109,106,106,109), MODOUT 

106 IF (OEL-DELMAXI 108,108,107 

107 DEL * M0DFIUEL,0ELMAX) 

108 IF (DEL) 114,109,114 

114 DELT * M IN1F ( DEL T , DEL ) 

109 IF ( XA8SF( lM00E)-4) 1,110,1 

110 CALL TUOES 

1M0DE * XSI0NFI2, IMOOE) 

1 CALL N80DY 

2 CALL EXTRAS 
C 

C PART 9. COMES HERE FOR END OF SU8 TRAJECTORY. 

IF (DONE) 113,111,111 

111 DONE * 0. 

IF (NSTAGE-LSTAGE) 112,115,115 

112 NSTAGE = NSTAG6+1 
GO TO 100 

113 DONE * 0. 

115 CALL EXTRA 

IF (RETURN) 103,116,100 

116 RETURN 
END 

SUBROUTINE STEP 

C SUBROUTINE STtP TESTS FOR THE END OF THE PROBLEM, COMPUTES STEP SIZE, AND 

C CONTROLS QUANTITY OF OUTPUT DATA. ENO OF PROBLEM OCCURS IF TIME = TMAX, 

C ST EPGQ+STEPnQ = STEPMX, OR C(LOOKX) * XLOOK. THE LAST OPTION ALLOWS STGP- 

C PING ON A DLPENUENT VARIABLE. THE TEST FOR STOPPING AT XLOOK IS NOT MADE 

C UNTIL C(LOUKSW) IS GREATER THAN SWLOOK. CONTROL UN QUANTITY OF OUTPUT IS 

C 

C M0D0UT=1 

C 

C 2 

C 3 

C 

C 4 

C 

COMMON C 
C 

DIMENSION A ( 600 ) , 8(700), C(4000), 

1 XPR I M ( 200 ) , OELTl (10) 

C 

EQUIVALENCE 

I ( A ,C ( 1 1 ) ) , ( A 1 ,8 ( 10)), (A2 ,8 

2 1 B ,C (1111)), (DEL MAX, A ( 19)), (DEL ,A 

3 ( DELT , B ( 1 ) ) , ( OONE ,B ( 39)), (E2 ,B 

4 (END , A ( 5)), (ERLOG ,B I 17)), (H2 ,B 

5 ( I NLOOK , A ( 599) > , (LOOKSW.A ( 9)),(L00KX ,A 

6 ( MQDOUT , A ( 20 ) ) , ( NSTAGE , A ( 3)), (OELTl ,A 

7 ( RAT 10 ,B ( 58) ), (SIGNAL, B ( 3 1 )),( SPACES , B 

8 ( STEPGO, A ( 41) ), (STEPMX, A ( 16) ) , ( STEPNO, A 

9 ( STEPS , A < 17) ),( SWLOOK, A ( 10)), (TABLE ,C 

EQUIVALENCE 

1 I TMAX , B ( 4)), (TWIN ,A ( 18)), (TTOL ,A 

2 ( XLOOK , A ( 12) ) , ( XPR IM ,C ( 711)),(XT0L ,A 

3 ( NSTART , B ( 24 ) ) , { SWI TCH, A < 60 1 > ) , < OUTPOT , B ( 

CHECKFt A,B,C) = ABSF(A-B) - ABSF(A-C) 

C 

C PART 1. TEST FOR ENO OF THE PROBLEM (MAXIMUM PROBLEM TIME 

C NUMBER OF STEPS). 

STEPGO = STEPGO ♦ 1. 

OUT * OUTPOT 

IF ( ABSF(TMAX-XPRIMIl) J-TTOL) 1,1,3 

1 OONE « 1.0 
112 CALL OUTPUT 

IF (OUTPOT) 26,111,26 
111 WRITE OUTPUT TAPE 6, 2, NSTAGE 

2 FORMAT ( 6 HOST AGE 1 2 , 1 1H COMPLETED.//) 

GO TO 26 

3 IF ( STEPGQ+STEPNO— STEPMX ) 7,4,4 

4 CALL OUTPUT 

WRITE OUTPUT TAPE 6, 5, STEPMX 

5 FORMAT ( 22H0STEPGQ*STEPN0=STEPMX*F6. > 

CALL EXIT 

C 



OUTPUT EVERY NTH STEP ( N*ST6PS ) UNTIL TIME = TMIN, THEN 
GO TO MO0E 2 . 

OUTPUT AT INTERVALS OF DELMAX UNTIL TIME = TMAX. 

OUTPUT AT INTERVALS OF DELMAX UNTIL TIME = TMIN, THtN 
GO TO MODE 4 . 

OUTPUT EVERY NTH STEP UNTIL TIME * TMAX. 


63 


C PART 2. COMPUTE STEP SIZE (OELT) AND CONTROL OUTPUT. 

7 N- 1 

A3 = ( A2-AI) -RATIO+A2 
AA = (ERL0G-A3I/5. 

IF { (ABSF(AA)-88.028) *ABSF(SWITCH) > 8,8,60 

8 OELT = SIGNF(£XPF( AA) , DbLT > 

IF (DELT/H2-3.) 10,10,9 

9 DELT = 3.*H2 

10 MODOUT = MOUOUT 

GO TO ( 11, li> , 13, 21 ) .MODOUT 

11 IF(DELT*(XPR1M( 1 ) + 3. *DEL T-TM I N ) ) 21,12,12 

12 MODOUT = 2 

DEL = TM IN — XPRIM( 1 > 

GO TO 16 

13 I F ( DEL T • (XPRIM(l) - THIN)) 15,15,14 

14 MOOOUT = 4 
GO TO 21 

15 DEL = DEL— HZ 

16 SPACES = IN IF ( !DEL/DELT)+SIGNF( .9, (DEL/DELT) ) ) 

17 I F { SPACES ) 20, 18,20 

18 CALL OUTPUT 
N=2 

DEL = DELMAX 

IF ( A8SF ( DEL ) - ABSF ( DEL T ) ) 19,16,16 

19 DELT = S IGNF ( DEL , OELT ) 

GO TO 16 

20 DELT = DEL/SPACES 
GO TO 23 

21 IF < M0DF1 STEPGO, STEPS) ) 23,22,23 

22 CALL OUTPUT 
N = 2 

C 

C PART. 3. SEARCH FOR CILOQKX) = XLOOK UNLESS L00KX=0. 

23 I F ( LOOK X) 27,42,27 

27 LOOK X = LOUK X 
LOOK SW = LOOK SW 
OUTPOT = 1. 

GO TO (44,45) , N 

44 CALL OUTPUT 

45 IF(SWITCH) 32,28,33 

28 I F ( SW LOOK - CILOOK SW)) 29,29,42 

29 XT0L1 = XTQL*ABSF( XLOOK) 

IF (XT0L1) 31,30,31 

30 XT0L1 = XTOL 

31 SWITCH = -1. 

GO TO 41 

32 SWITCH = 1. 

ASSIGN 43 TU MODE 
OVER = 0. 

F = 0. 

T=0. 

33 SLOPE * (CIL00KX)-0LDX)/H2 
GO TO MODE, (43,35) 

43 I F ( SLOPE * ( C 1 LOOK X) - X LOOK)) 350,41,41 
350 ASSIGN 35 TO MODE 

35 I F ( ABSF(C( LOOK X)- X LOOK) - XT0L1) 36,36,37 

60 T - 1 - 

36 IF (OUT) 63,46,63 

46 OUTPOT = 0. 

CALL OUTPUT 

63 IF (T) 61,47,61 

61 IF (OUT) 62,51,62 

51 WRITE OUTPUT TAPE 6,64, LOOK X , C ( LOOK X >, H2 , LOOK X , SL OP E 

64 FORMATt 3H0CI 14, 4H) = 1PG15.8.31H CONVERGENCE TROUBLE. DELT 
1G15.8, 14H SLOPE OF C(I4,13H) VS. TIME = G15.8//J 

GO TO 62 

47 IF (OUT) 62,50,62 

50 WRITE OUTPUT TAPE 6, 48, LOOK X, C(LOOK X) 

48 F0RMAT13H0CI 14 , 2H ) = 1 PG1 5. 8// ) 

62 LOOKX = 0 
XTOLl = 0. 

SIGNAL = 1. 

SWITCH = 0. 

DONE = END 
NSTART = 0 
NSTAGE=NSTAGE 

DELT = DELT1 ( NSTAGE ) 

49 CALL INPUT ( INLOOK, C, TABLE) 

IF ( OONE ) 110,42,110 

110 IF (OUT) 26,111,26 

37 SIGN = CHECKFIOLOX, XLOOK, C(LOOK X)) 

IF(SIGN) 40,40,38 

40 OVER = 1. 

GO TO 400 

38 IF (OVER) 400,401,400 

401 XGUESS * CILOOKX )+SLOPE*DELT 

IF (CHECKF(C(LOOKX), XLOOK, XGUESS)) 402,41,41 

402 F = F+l. 

IF IF-7.) 400,400,403 

403 SLOPE = SLOPE/F 

400 IF (SLOPE) 404,60,404 

404 DELT = SIGNF(ABSF( XLOOK-C (LOOKX) )/SL0PE,SIGN*H2) 

41 OLDX = C ( LOOK X) 

42 IF (ABSF(TMAX-XPRIMI 1) ) -AB SF ( DEL T ) ) 25,26,26 

25 DELT = T MA X— XPR I M ( 1) 

GO TO ( 26,24,24,26) , MODOUT 

24 DEL = DEL-DELT 

26 OUTPOT = OUT 
RETURN 

END 
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c 

c 

c 

c 

c 


c 


c 

c 


c 

c 


c 

c 


c 

c 


c 


SUBROUTINE STOAT A 

THIS ROUTINE CLEARS THE A, XPRIM, XPRIMB ARRARYS AND LOADS A SET DF 
STANDARD DATA INTO THE MACHINE. ANY VALUES SET HERE MAY BE OVERWRITTEN BY 
INPUT I IN THE MAIN PROGRAM. 

COMMON C 


DIMENSION AI600), B1700), C(4000), 



1 PNAME 112), 


AMASS (30), 

XPRIM (200), 


2 


COEFN (190), 

ICC 14), 


3 AK (31, 


XDOT (100), 

IND (3), 


4 REFER (12) 

5 RMASS1 (10) 

* 

RCR1T (30), 

AW (4), 


equivalence 

HA ,C ( 

11) ) 

,<AK ,A ( 

51)), (AMASS , A 

( 347)) 

2 ( AU , A I 

29) ) 

, ( AW , A ( 

55)), IB ,C 

(1111) ) 

3 ( BODYCD , A ( 

143} ) 

.(COEFN , A { 

407) ) , (CONSTU,A 

( 32)) 

4 ( CONSU , A ( 

31) ) 

t ( DTOFF J , A ( 

2 3 ) ) , ( EREF , A 

( 13) ) 

5 ( ERL I MT , A ( 

14) ) 

t I E TOL * A ( 

30) ) , ( GASFAC , A 

( 46) ) 

6 ( I CC , A ( 

153) ) 

, ( IMOOE , A ( 

1 ) ) t ( I ND • A 

( 60) ) 

7 ( LOOKSW , A ( 

9) ) 

, IMUDOUT.A ( 

20 ) ) , ( NEQ , A 

( 2) ) 

8 ( N S T AGE ? A ( 

3) ) 

, ( OBLATD, A ( 

27) ) , (OBLATH, A 

{ 28) ) 

9(08LATJ,A ( 

26) ) 

.(PNAME , A ( 

287) ) , (RCR1T , A 

( 377)) 

EQUIVALENCE 
KREFER , A ( 

317) ) 

• (RE * A ( 

25)),(SPD , A 

( 44) ) 

2 ( SQRDK 1 , A ( 

47) ) 

, ( STEPMX , A ( 

16) ),( STEPS , A 

t 17) ) 

3( TFILE * A ( 

6) ) 

•(XDOT , 8 ( 

501) ) , (XPRIM ,C 

( 711)) 

4 ( XTOL , A < 

11) ) 

, ( RMASS1 • A ( 

73) ) 



CLEAR INITIAL CONDITIONS AND CONTROL PARAMETERS. 
DO I J=l,1100 
I A { J ) = 0. 


THE FOLLOWING NH STATEMENTS LOAD THE BODY NAMES INTO THE 

PNAHEt 1 ) 

= 

3HSUN 

PNAME12} 

= 

6HMERCUR 

PNAME ( 3 ) 


5HVENUS 

PNAMEI 4 ) 

= 

5HEARTH 

PNAME ( 5) 

= 

4HMARS 

PNAME ( 6 ) 


6HJUPITE 

PNAMEI 7) 

= 

6HSATURN 

PNAME ( 8 ) 

= 

6HURANUS 

PNAME ( 9 ) 

= 

6HNEPTUN 

PNAME ( 10 1 

■ 

5HPLUT0 

PNAHEt ID 

* 

4HM00N 

PN AM6 ( 12 ) 

3 

6HEARTHM 


MACHINE. 


FILL OUT SUN REFERENCE LIST. INITIALIZE MASS ARRAY. 
DO 2 K=1 , 10 
RMASSUK) * 1. 

2 REFERIK+l) = PNAMEI l ) 

REFER ( 12 ) = PNAMEI l ) 


FILL OUT EARTH REFERENCE LIST. 
REFER! 1 ) = PNAME ( 4 ) 

REFER! 4) * 5HZER0+ 

REFERl 11 ) * PNAMEI4) 
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C LOAD THE REMAINING STANDARD DATA. 

AMI) = 0.5 
AKI2) * 0.5 
AK ( 3 ) = 1.0 
AMASS! 1) = 1.0 
AMASS! 2) = I. 0/6120000.0 
AMASS! 3) = 1.0/408645.0 
AMASS! 4) = 1.0/332951.3 
AMASS! 5 ) = 1.0/3088000.0 
AMASS ( 6 ) = 1.0/1047.39 
AMASS ( 7 ) = 1.0/3500.0 
AMASS ( 8 ) = 1.0/22869.0 
AMASS { 9 ) = 1.0/18889.0 
AMASS! 10) = 1.0/400000.0 
AMASS (11) =AMASS(4)/81.335 
AMASS! 12) =AMASS ( 4 ) + AMASSlll) 

AU = 1.49599 Ell 
AW! 1 )=l./6. 

AW(2)=AW! 1)+AW!1) 

A W 1 4 ) = AH ! 1) 

AW(3)=l.-( Aw!2)+(AW( 1) +AW14) ) ) 
80DYCD = PNAME ( 4 ) 

COEFNf 1) = -1E20 
COEFN! 189) = 1E20 
CONSTU = 1.0 E-6 
CONSU = IE— 6 
ETOL = 0.01 
DTOFFJ = 244. E4 
EREF = IE— 6 
ERLIMT = 3E-6 
GASFAC * 20.064881 
ICC! 1 ) =185 
ICC ( 2 ) =185 
ICC ( 3 ) =185 
ICC! 4) =185 
IMODE = 1 
IND! 1)=2 
IND! 2)=3 
IND! 3 ) =1 
LOOKS W = 711 
MOOOUT = 4 
NEQ=8 

NSTAGE = 1 

08LATJ = 1.62345 E-3 
O0LATH = -5.75 E-6 
OBLATD = 7.875 E-6 
RCRIT(l) = 1.0 E+20 
RCRITI2) = 1.0 E + 8 
RCR I T ( 3 ) = 6.14 E+8 
RCRITI4) = 9.25 E+8 
RCR l T ( 5 ) = 5.78 E+8 
RCR I T ( 6 ) = 4.81 E+10 
RCR I T ( 7 ) = 5.46 E+10 
RCR IT I 8 ) = 5.17 E+10 
RCR I T ( 9 ) = 8.61 E+10 
RCR IT! 10) =3.81 E+10 
RCR I T ( 1 1 ) =1.60 E+8 
RE = 6378165. 

SPD = 86400.0 

SQRDK1 = 2.959122083 E-4 

STEPMX= 100. 0 

STEPS = 1. 

TFILE = 1. 

XDOT(l) = 1.0 

XPR I M ( 2 ) = RMASSl(l) 

XTOL = 5E-8 
WRITE OUTPUT TAPE 6,3 
3 FORMAT ! 15H0ST ANDARD DATA.) 

RETURN 

END 
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SUBROUTINE TESTTR 

C SUBROUTINE TESTTR MAY BE CALLED FOR ONE OF TWO REASONS, 111 TO TEST FOR AND 

C POSSIBLY TRANSLATE THE ORIGIN (WHEN IMODE IS ♦ ) OR (2) TO CHANGE THE 

C VARIABLES OF INTEGRATION (WHEN IMODE IS A TRANSLATION OF THE ORIGIN 

C OCCURS WHEN THE OBJECT MOVES INTO A SPHERE OF INFLUENCE WHICH IS SMALLER 

C THAN ANY OTHERS IT MAY ALSO BE IN. WHEN THIS HAPPENS, THE NAME OF THE NEW 

C ORIGIN IS MOVED TO THE BEGINNING OF THE BNAME LIST AND ORDER IS 

C CALLED TO REORDER THE BNAME LIST. 

C 

COMMON C 
C 

DIMENSION A ( 600 ) , BI700), CI4000), 

1 XPRIM( 100,2) »XPRIMB( 100,2) ,XWH0LE(6) , VEFMI 3,8 ) , VXI 3) , 

20R8ELS ( 6) , BMASS ( 8),BNAME(8),RB(3,6), RBCRITt 8 ) , R( 8) 


EQUIVALENCE 


1 ( A ,C 

( 1 1 ) ) , ( AMC , B 

( 87) ), (ASYMPT, A 

( 7)), 

2 ( B ,C 

(lllU ), (BMASS ,B 

( 137) ) , (BNAME ,B 

( 122)), 

3 ( CHAMP ,B 

I 25 ) ) , I DELT ,B 

( 1)) , ( GK2M , B 

( 36)), 

4( IMODE , A 

( 1) ), (NBODYS, B 

( 41) ), (ORBELS, B 

( 116)), 

5 ( RBCR I T , B 

( 145)), (RB , B 

l 193 ) ) , ( REVS , A 

£ 48)), 

6 ( R , B 

( 102)), (SQRDK ,B 

I 35)), (TABLE ,C 

(1911) ) , 

7ITMAX ,B 

( 4) ), (TRSFER, B 

( 8)), (TRU ,3 

( 40)), 

B l TTEST , A 

( 54 ) ) , ( VEFM ,B 

( 241)), (VX ,B 

( 92)), 

9IXPRIM ,C 
EQUIVALENCE 
1 ( OUTPOT , B 

( 711) ) , (XPRIMB,C 
( 399)) 

( 911) ), (XWHOLE, B 

( HO)) 


IMODE = IMODE 
IF (IMODE) 12,12,1 
C 

C IF IMODE IS +, TEST FOR TRANSLATION OF THE ORIGIN. 

1 CHAMP= l.E+30 

DO 4 JB= L » NBODYS 
IF (R(JB)— RUCRIT(JB) ) 2,4,4 

2 IF ( CHAMP— R BCR IT ( JB ) ) 4,4,3 

3 CHAMP = RBCRITIJB) 

NCHAMP = JB 

4 CONTINUE 

IF ( NCHAMP— 1 ) 26,26,5 

5 TRSFER = 1.0 

8 BTEMP = BNAME ( 1 ) 

BNAME ( 1 ) = BNAME ( NCHAMP ) 

BNAME ( NCHAMP ) = BTEMP 
TTEST = 0. 

REVS = 0. 

IF (OUTPOT) 6,9,6 

9 WRITE OUTPUT TAPE 6 , 10 , BNAME ( NCHAMP ), BNAME ( 1 ) 

10 FORMAT ( 28H00R IG IN IS TRANSLATING FROM A6,4H TO A6) 

6 CALL EPHMRS 
DO II K= 1 , 3 

VX(K> - VX(K)-VEFMIK, NCHAMP) 

RB ( K ) * RB(K, NCHAMP) 

XPR1M(K+2,1)=VX(K) 

XPRIM(K+5* 1 ) =R B ( K ) 

XPR I MB ( K+2 , 1 ) =* 0. 

XPR I MB ( K + 5 » 1 ) - 0. 

XWHOLE ( K ) = VX( K ) 

11 XWHOLE ( K + 3 ) = RB ( K ) 

GO TO 20 

C 

C IF IMODE IS -, CHANGE THE VARIABLES OF INTEGRATION. 

12 DO 13 K-1,3 
XPRIMtK+2, 1) =XWHOLE( K) 

XPRIMI K + 5 » 1)=XWH0LE( K+3) 

XPRIMBf K+2* 1 ) * 0. 

XPR I MB ( K + 5 , 1 ) a 0. 

VX(K) = XWHOLE ( K ) 

13 RBIK) = XWHuLE { K+3 ) 

GO TO (16,14,15) , IMODE 

14 COOE = 5H0RBIT 
IMODE = 1 

GO TO 18 

15 IMODE * 3 
GO TO 17 

16 IMODE *= 2 

17 COOE = 6HRECTAN 

18 NCHAMP = 1 

IF (OUTPOT) 20,7,20 

7 WRITE OUTPUT TAPE 6, 19, CODE 

19 FORMAT ( 33HO INTEGRAT ION MODE IS CHANGING TO A6 ) 

20 GO TO ( 21,26,26) , IMODE 

21 CALL C0NVT1I VX,AMC) 

GK2M= SQRDK* ( BMASS ( NCHAMP } +XPR IM ( 2, 1) /1.9866 E + 30) 

30 CALL C0NVT2 

C IF ORIGIN TRANSLATION CAUSES PATH TO LIE NEAR AN ASYMPTOTE, CHANGE 

C INTEGRATION VARIABLES TO RECTANGULAR IF THEY ARE ORBIT ELEMENTS. 

IF ( QRBELS ( 1 )— 1 . ) 24,24,22 

22 IF ( ABSF ( TRU)— 2. 3/ SQRTF ( ORBELS ( 1 ) ) ) 24,24,23 

23 ASYMPT = 1.0 
GO TO 15 

24 DO 25 J= l , 6 

25 XPRIMIJ+2,1) * ORBELS ( J) 

26 IF (TRSFER) 27,28,27 

27 CALL INPUT ( 101 , C, TABLE ) 

29 CALL ORDER 

28 RETURN 
END 



SUBROUTINE THRUST 

C THIS ROUTlNt COMPUTES X,Y,Z THRUST ACCELERATIONS. THE THRUST VECTOR IS 

C ASSUMED COINCIDENT WITH THE LONG I TUND I NAL AXIS OF THE VEHICLE, WHICH IS 

C ORIENTED TO THE RELATIVE WINO VELOCITY BY THE ANGLE OF ATTACK (ALPHA) AND 

C THE ROLL ANGLE (BETA). ALPHA IS ASSUMED TO BE A QUADRATIC FUNCTION OF TIME 

C WHEREAS BETA IS ASSUMED TO BE CONSTANT. 

C RE VOL V IS THE EARTHS ROTATIUN RATE IN RADIANS/SEC ( 7. 292 l 1 585E-5 ) AND THE 

C FACTOR 8589934592.= 2**33 IS REMOVED TO PREVENT OVERFLOW. 

C 

COMMON C 


C 


c 


c 


DIMENSION A ( 600 ) « B(700), CI4000), 

I FORCE ( 3 ) , PAR ( 3 ) , VATM(3), P(3), IND( 3 ) , RAMC ( 5 ) , RB( 3 ) , X ( 100 ) 


EQUIVALENCE 

1 ( A ,C 

2 ( B ,C 

3 ( COSBET , B 
4 ( FORCE ,B 
5 ( PMAGN ,B 
6IPUSH0 ,B 
7 I R ATMOS , B 
8 ( R ,B 

9 ( S INALF, B 

EQUIVALENCE 
l( VQ , B 
2 ( V Y , B 


( 11 ) ) , t AEXIT ,B 
( 1111) ), (BETA , A 
( 49) ) , (EXITA ,8 
( 66)), I I NO , A 
( 50)), (PRESS ,B 
( 391)), (PUSH , A 
( 23)), (R8 , B 
( 102) ) , ( RSQRD ,B 
( 46) ) , ( SINBE T , B 

( 100) ) , ( VQSQRD,B 
l 93)), (VZ , B 


( 3)), (ALPHA , A 

( 50 ) ) , ( CQSALF, B 

( 392)), (FLOW , B 
( 60)), (PAR ,B 

( 33)), (P , B 

( 166)), (RAMC , B 
( 193) ) , (REVOLV.B 

( 45)), (SIMP , B 

( 47) ) , ( VATM ,B 

( 10 L ) ) , ( VX , B 
( 94)), (X , B 


49) ), 
48) ) , 
5)), 
60) ), 
84) ), 
393) ) , 
21 ) ) , 
2 ) ), 
97) ) 

92) ) , 
401) ) 


SINBET = SINF( BETA/57.2957795) 

COSBET = COSFIBETA/57. 2957795) 

VATH ( l)=VX+REV0LV*RB(2) 

VATM( 2)=VY— REVQLV*RB( 1 ) 

VATM( 3)=VZ 

3 CALL CONVTll VATM, RAMC) 

4 ALPHA = QUAO(X( 1 ), 1 )/57. 2957795 
SINALF=SINF( ALPHA) 

COSALF=CQSF( ALPHA) 

DO 1 J 1- 1 » 3 
J2 = I NO ( J 1 ) 

J3= I ND ( J2 ) 

1 P(J1) * (VATM(J2)* RAMC (J3)-VATM(J3) *RAMC ( J2 ) ) / 8589934592 . 

PMAGN= SQRTF ( P ( 1 )*P(1)+P(2)*P(2)+P(3)*P(3)> 

PUSH = PUS HO— EX ITA*PRESS 
TDPMAG = PUSH/ PM AGN/X 1 2 ) 

R4 = SINBET/VQ 

R5 = C0SALF/RAMC14) 

DO 2 J 1=1, 3 
J2= IND ( J 1 ) 

J3= I ND ( J2 ) 

PARI J1)=P( J2)*VATM( J3)-P( J3)*VATM( J2) 

2 FORCE ( Jl ) = TDPMAG*( SINALF*(COSBET*P(Jl) +R4*PAR ( J 1 ) ) — R5* ( P ( J2 ) * 

1 RAMCt J3)-P( J 3 ) *RAMC( J2 ) ) ) 

RETURN 

END 


SUBROUTINE TUDES 

C THIS ROUTINE COMPUTES THE RECTANGULAR POSITION AND VELOCITY COMPONENTS 

C WITH RESPECT TO THE EARTH MEAN EQUINOX AND EQUATOR OF 1950.0 FROM THE 

C LATITUDE, LONGITUDE, AZIMUTH, ELEVATION, ALTITUDE, TOTAL VELOCITY, AND 

C TIME. ALSO, WHEN TKICK DOES NOT EQUAL ZERD, A NON-DRAG VERTICAL STEP OF 

C SIZE TKICK IS MADE IN CLOSED FORM (STATEMENTS 2 TO 4). THE INTEGRATION 

C WILL THEN BEGIN AT TIME EQUAL TO TIM6+TKICK WITH THE ORIENTATION SPECIFIED 

C BY THE ABOVt FOUR ANGLES AND THE COMPUTED VALUES OF ALTITUDE AND VELOCITY. 

C FDR THE CLOSED FORM APPROXIMATION, A CONSTANT FLOW RATE (FLOW), VACUUM 

C SPECIFIC IMPULSE (SIMP) AND ENGINE EXIT AREA (AEXIT) ARE ASSUMED KNOWN. 

C THE ATMOSPHERIC PRESSURE IS TAKEN TO BE THE SEA LEVEL VALUE. 

C 

COMMON C 
C 

DIMENSION A l 600 ) , B(700), CI4000), 

1 S I N A ( 4 ) , COSA ( 4 J , ANGLEBI 4 ) , XPRIM(200) 

C 


EQUIVALENCE 


I ( A ,C 

( ID), (AEXIT 

, B 

( 3)), (ALT , A 

( 4) ) 

2 ( AZ I , A 

( 35)), (B 

,c 

(1111) ) , ( DTOFF J , A 

( 23) ) 

3 ( EL E V , A 

( 36)), (FLOW 

, B 

( 5 ) J , ( GK2M , B 

( 36) ) 

4 ( L AT , A 

( 33)), (LONG 

, A 

( 34) ) , ( OBL ATJ , A 

( 26) ) 

5 ( OBL ATN » A 

( 40)), (RE 

, A 

( 25) ) , (RESQRD.B 

( 7) ) 

6(R0TATE, A 

( 39)), (SIMP 

, B 

( 2) ), (SPD ,A 

( 44) ) 

7 ( STEPGO, A 

( 4 1 ) ) » ( S TEPNO 

, A 

( 42) ) , (TKICK , A 

( 15) ) 

8 ( TOFFT , A 

l 24)),(VEL 

,A 

( 37 ) ) , ( XPR I M ,C 

( 711)) 

9 ( OUTPOT , B 

1 399)) 





EQUIVALENCE (QLAT , LAT ) , ( QLQNG , LONG ) 


C 
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ALT1 * 0. 

VEL1 » VEL 
DELI * 0. 

DEL = 0, 

ASSIGN 1 TO NGO 

OAYS = DTOFFJ - 2433282.5 

GREEN = MOUFI 100. 0755426+. 98 5647346DAYS+2 . 90 1 5E- 1 3DA YS**2 
1+7.29211585E— 5*( TOFFT* SPU+XPR IM ( 1) 1*57.2957795,360.) 

SINA(i) * SINFtQLAT/57. 2957795) 

IF (OBLATN) 102,101,102 

101 RADIUS * RE + ALT 
GO TO 8 

102 RADIUS=6356783.20/SQRTFI .9933065783+. 00669342 1 685*S INA I 11**21 + ALT 
GO TO 8 

1 XPRIHI6) * C0SA{2)*C0SA( 1 1 * RAD I US 
XPR I M( 7 1 = SINA(2) *COSA( 1) *RADIUS 
XPR I M ( 8 1 = SINA( L)*RAOIUS 
RMASSO = XPRIMI2) 

XPR I HI 2 1 - XPR I H ( 2 )— FLO W*TK ICK 
IF (OUTPOT) 12,11,12 

11 WRITE OUTPUT TAPE 6 , 3, STEPGO, STEPNO, LAT , LONG, A2 1 , ELEV, ALT , XPR I MI 
11) , VEL, RMASSO, ( XPR IM t J ) , J=6 , 8 ) 

3 FORMAT ( 6H0STEP=F5. *2H +F4.,4X,6H LAT. = 1PGL 5. 8 , 7H L0NG.=G15.8 , 6H AZ 
1 1 • =G15. 8 , 7H ELEV. = G15.8,6H ALT. = G1 5. 0/6H T IME = Gl 5 . 8 , 6H VEL.=Gl5-8, 
67H RMASS=G15.8,4X,2HX=G15.8, 5X , 2HY=G15 . 8 , 4X , 2HZ=G1 5. 8 ) 

12 IF (TKICK) 2,50,2 

2 XPRIM(l) = XPRIMt 1 J+TKICK 
B1 = LOGF I RMASSO/ XPR IM I 2 ) ) 

SIMPSL = S IMP— AEX I T/ FLOW* 10332.275 
VEL 1 = VEL+S IMPSL+9. 80665*81— G*TK ICK 

ALT1 = TK 1 CK* ( VEL— G*TK ICK/ 2. +9. 80665* S IMP SL* ( 1.-B1+XPRIM12)/ 
l I RMASSO— XPR IM I 21)1) 

4 RADIUS = RAUIUS + ALT 1 

GREEN = GREEN + 7. 2921 1585E-5*TK I CK*57. 2957795 
ASSIGN 5 TO NGO 
GO TO 8 

5 XPR I MI 6 ) = G0SA(2)*C0SA{ 1|*RAUIUS 
XPR I MI 7 1 - SINA(2)*C0SA( l) -RADIUS 
XPRIMIB) a SINAI 1) -RADIUS 

50 IF I OBLATN) 6,7,6 

6 DELI = ATANF I ( C2— 1. ) / I C3— 1. I*SINA( ll/COSAI 1) 1*57.2957795— QLAT 

7 DEL2 = RAD1US/G*SINAI 1)*C0SA( 1 1 *ROTATE*ROT ATE* 57 . 295 7795 L 
DEL = DELI + DEL 2 

ASSIGN 10 TO NGO 

8 ANGLEBU) = QLAT + DEL 
ANGLEBI2) = QLONG + GREEN 
ANGLES 1 3) = AZI 

ANGLEBI 4 ) * ELEV 
DO 9 1*1,4 

SINA(I) = SINFI ANGLEBI I 1/57.2957795) 

9 COSA(I) = COSFIANGLEBI I)/57. 2957795) 

Cl = 5.*RESQRD/RADIUS/RADIUS*0BLATJ 
C2 = Cl* (SINAI 1 ) *SINAI 1 ) — .6) 

C3 = C 1* I S INAI 1 ) • SINAI 1 )— . 2 ) 

G = GK2M/RAU IUS/RAD1US 
GO TO NGO, (1,5,10) 

10 COSl - COSAl 1)*SINAI4)-C0SAI4)*C0SA(3)*SINA( 11 
COS 2 = COSAl 4) *SINA( 31 

XPR I M I 3 ) = VEL1*ICOS1*COSA(2)-COS2*SINA(2) 1 -XPR I M { 7 1 *R0T ATE 
XPR IM ( 4 1 = VEL1*(C0S1*SINA(2) +C0S2*C0SA (2) )+XPRIM(6)*R0TATE 
XPRIMI5) = VEL1*ISINAI 1)*SINA(41+C0SA( 1)*C0SA(3)*C0SA(4) ) 

RETURN 

END 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE TAPE 

SUBROUTINE IAPE USES THE MASTER MERGED E PHEMER I DES TAPE (TAPE 9 AT LEWIS) 
TO COMPILE A WORKING EPHEMERIS TAPE (TAPE 3 AT LEWIS) WHICH CONTAINS ONLY 
THAT OATA NtEOED AT EXECUTION TIME. THIS MINIMIZES TAPE HANDLING DURING 
EXECUTION- 2 EPHEMERIS FILES ARE ON TAPE 9, FIRST FILE HAS DATA AND IS 
IDENTIFIED BY THE SECOND WORD OF EACH 254 WORD RECORD (FIRST WORD IS THE 
DUMMY FORTRAN COMPATIBLE WORO, SECOND W0RD=2 1 . THE SECOND FILE IS ONLY 2 
WORDS LONG, FIRST WORD IS FORTRAN COMPATIBLE, SECOND W0RD=3). 

MASTER FILE 1 — PLANETS (EXCEPT MERCURY AND EARTH), SUN, MOON, AND 

EARTH-MOON BARYCENTER FROM SEPT. 25, 1960 TO ABOUT 2000. 
EACH EPHEMERIS COMPILED REQUIRES A SET OF INPUT 300 DATA. THE FIRST PIECE 
OF DATA WRITTEN ON A FILE IS THE FILE IDENTIFICATION NUMBER, FILE. EACH 
FILE IS NUMBERED CONSECUTIVELY STARTING WITH FlLE=l. SINCE MUUN DATA IS IN 
TERMS OF EARTH RADII, THE CONVERSION OF MOON DATA TO A.U. IS MADE BEFORE 
WRITING ON TAPE 3. THE COMMON USED IN SUBROUTINE TAPE IS LOCAL AND ALL 
BUT TAPE3 IS CLEARED BY A FINAL CLEARING LOOP. 

FUNCTION COMPARF I A , B ) IS EQUIVALENT TO (A-B) BUT WILL NOT OVERFLOW. 

NORMAL INPUT - EL I ST , TBEGIN, TEND, TAPE3 

EL I ST— THE BCD LIST OF EPHEMERIS DATA NAMES TO BE PLACED ON 

TAPE 3 . THE NAMES ARE READ FROM CAROS, AND IS USED TO 
MAKE THE TMAKE LIST. EL I ST IS NOT CHANGED IN STORAGE UNTIL 
THE FINAL CLEAR FOR THIS SUBROUTINE. 

TMAKE- THE LIST OF EPHEMERIS NAMES WITH DUPLICATES DROPPED AND 

ZERO SPACES CLOSED IN. AS THE EPHEMER I DES ARE FINISHbO THE 
NAMES ARE ERRASED FROM THIS LIST. 

TMADE— LIKE TMAKE BUT IS HELD FOR OUTPUT. 

TBEGIN- THE BEGINNING DATE EXPRESSED AS A JULIAN DAY. 

TENO- ENDING DATE EXPRESSED AS A JULIAN DAY. 

INTVAL— THE APPROX. NUMBER OF OAYS COVERED BY ONE SET OF COEFF. IT 
IS USED TO DECIDE WHICH DATA ARE TO BE ENTERED DOUBLE. THE 
DOUbLE ENTRIES PERMIT FASTER OPERATION IF REVERSAL OF 
INTEGRATION IS REQUIRED FOR ANY REASON. 

EDATE- JULIAN ENDING DATE FOR THE MASTER EPHEMERIS. 

ERTGAU— EARTH RADII PER A.U. 

COMMON C 
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DIMENSION 

1 C 1700)# TMAKE (12), LIST (30), 

2 EDATE (12), INTVAL (30), KTAG (12), 

3 EL 1ST 111), THADE (I2) r INTVA (2), 

A PNAME (30), T DATUM (252), DATUMT (21,12) 

C 

EQUIVALENCE 

1( TAPE3, C ( 2 ) ) , ( ERTOAU, C ( 3 1 > , I KTAG, C ( 4)),( FILE,C( 16)), 

2 ( EL IST»C( 17) ) , ( TBEGIN, C( 29)), ( TEND , C ( 30)), ( PNAME, C( 31)), 

3 ( KHAMP, C ( 61)), ( TMAD£,C( 73)), ( TMAKE, C( 85 > ) , ( TDATUM, C 1441 ) ) , 
4( EDATE, C( 127) ),( INTVAL, C( 157) ),( INTVA , C( 156 )),( DATUMT , C ( 1 89 ) ) 

C 

B COMPARE ( A, ti ) = ( A+B ) * ( — ( A*B) ) 

REMIND 3 
DO 1 K*1 » 4000 

1 C(K) = 0-0 
C 

C THE FOLLOWING NH STATEMENTS LOAD THE BODY NAMES INTO THE MACHINE. 

C NOTE. THE EARTH IS NOT IN THIS LIST (NO EPHEMER I S FOR EARTH.) 

PNAME ( 1 ) = 3HSUN 
PNAME ( 2 ) = 6HMERCUR 
PNAME ( 3 ) = 5HVENUS 
PNAME ( A ) = AHMARS 
PNAME ( 5 ) = 6HJUPITE 
PNAME ( 6 ) = 6HSATURN 
PNAME ( 7 ) = 6HURANUS 
PNAME ( 0 ) = 6HNEPTUN 
PNAME ( 9 ) = 5HPLUTU 
PNAME(10>= AHMODN 
PNAME ( 11)= 6HEARTHM 
C 

C PART 2- SET UP JULIAN DATES ENDING EACH EPHEMER I S. 


EDATE(l) = 2A51872- 5 11/2A/0Q 
EDATE( 3 ) = 2A518A8.5 10/31/00 
EDATE ( A ) = 2A51020- 5 7/26/90 
EDATE ( 5 ) - 2A73520- 5 2060 
EDATEI6) = 2A73520. 5 2060 
EDATE ( 7 ) = 2A73520- 5 2060 
EDATEt 8) = 2473520.5 2060 
EDATE ( 9 ) = 2A73520- 5 2060 
EDATE! 10 ) = 2AA0916- 5 11/26/70 
EDATEt 11)= 2A518A8.5 10/31/00 


INTVA = 30000 
INTVAL ( l) = 0 
INTVAL ( 2 ) * 5 
INT VAL ( 3 ) - 15 
INTVAL ( A ) = AA 
INTVAL ( 5 ) = 330 
INTVAL16) = 825 
INT VAL ( 7 ) = 1211 
INTVAL ( 8 ) = 1172 
INTVAL ( 9 ) = 1101 
INTVALt 10) = 2 
INTVAL (11) = 15 
FILE = 1. 

ERTOAU = A.265A6512 E-5 

2 MOON = 0 
LI = 1 

C 

C PART 2B- CALL INPUT AND SEE IF TAPE IS TO BE MADE- INPUT MUST ALWAYS 

C MAKE TAPE3=0.0 IF TAPE IS TO BE MADE. 

TAPE3 = 3. 

8 CALL INPUT (300, C, LIST) 

IF ( TAPE3 ) 63,3,63 

3 IF (FILE-1.) 20,10,20 

10 CALL SKFILE(9,2) 

C 

C PART 3. TAPE IS TO BE MADE SO MOVE EPHEMERIS LIST TO TMAKE ANO 

C TO TMADt (FOR OUTPUT), CANCEL ANY ZERO OR DUPLICATE NAMES. 

20 KOUNT = I 
DO 6 K= 1 » 1 1 
TMAKE ( K ) = 0. 

TMADE1K) = 0. 

A DO 5 J=l, KOUNT 

IF (COMPARF(ELISTIK) , TMAKE! J-l) ) ) 5,6,5 

5 CONTINUE 

TMAKE ( KOUNT ) = ELIST(K) 

TMADE ( KOUNT ) = ELIST(K) 

KOUNT = KOUNT+1 

6 CONTINUE 

KOUNT = KOUNT - 1 
C 

C PART 4. FIND INPUT ERRORS. 

7 IF(TB EG IN- 2437202-5) 66,9,9 

9 KM = 2 

11 ERROR «= 0. 

WRITE TAPE 3, FILE 
00 21 J=l, KOUNT 
KTAG(J) = 0 

12 DO 13 K- l ,20 

IF (COMPARF(PNAMEIK) , TMAKE! J) ) ) 13,16,13 

13 CONTINUE 


C PART 5. PRINTS OUT THE MISSPELLED NAMES AND OTHER ERRORS. 

14 PRINT 15, TMAKE ( J ) , TBEGIN, TEND 

WRITE OUTPUT TAPE 6 , 15, TMAKE ( J ) , T8EGIN, TEND, ( PNAME ( K ) , 
1EDATE(K) , K= l , 20) 

15 FORMAT! 23H TROUBLE ON TAPE 3 MAKE / 2X,A6,10H T 8EGIN= F10.1.BH 
l T END= F10.1//2(2X,A6,F20.1) ) 

ERROR = 1. 

GO TO 21 
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PART 48. CHtCKS DATES AND STORES INDEX FOR MOON SO THAT EARTH 
RAO 1 1 CAN BE CONVERTED TO A.U. 

16 IF f 10— KJ 10,17,18 

17 MOON = J 

10 KTAG(J) * K 

19 IF (EOATE(K)- TEND) 14,21,21 

21 CONTINUE 

ASSIGN 36 TO NS1 
IF (ERROR) 22,2 2,68 

PART 6. FIX UP A TAG ( KTAG) TO INDICATE WHETHER TO ENTER DATA DOUBLE OR 
NOT. KHAMP WILL BE SHORTEST INTERVAL. KTAG WILL BE NON-ZERO IF 
ANY DATA ENTERS MORE THAN ONCE FOR 10 ENTRIES OF THE MOST 
FREQUENT DATA. 

22 KHAMP * INTVAL(O) 

DO 23 J=l, KOUNT 
K - KTAG(J) 

KHAMP » XM INOF (KHAMP, INTVAL(K) ) 

23 CONTINUE 

KHAMP = KHAMP *10 
DO 24 J= 1 , KOUNT 
K * KTAG(J) 

24 KTAGU) - INTVAL(K) / KHAMP 

PART 7. LOCATE FILE 2 ON TAPE 9. 

25 READ TAPE 9, KFILE 

26 IF (KM-KFILE) 27,31,29 

27 IF (KFILE - 3) 20,20,29 

28 BACKSPACE 9 
BACKSPACE 9 
CALL BSF ILE( 9) 

GO TO 25 

BY PASS A FILE. 

29 CALL SKF I LE ( 9 ) 

GO TO 25 

PART 0. THIi IS CORRECT FILE ON TAPE 9, READ DATA. THERE CAN BE UP 
TO 12 SETS OF DATA PER RECORD. A SET OF DATA IS 21 WORDS. 

31 BACKSPACE 9 

32 READ TAPE 9, KTAPE, ( TOATUM(I) , 1*1,252) 

GO TO NS 1 , (36,46) 

PART 9. IS THIS A SATISFACTORY STARTING POINT, QUESTION MARK. 

THE 1ST SET OF DATA FOR EACH PLANET MUST PRE DATE TBEGIN. 

PART 9 IS EXECUTED ONLY ONCE. 

36 DO 42 J=L I , KOUNT 
DO 37 K* 1 , 232 , 2 1 

IF ( COMPARE (TDATUM(K), TMAKE ( J ) ) ) 37,39,37 

37 CONTINUE 

38 LI = J 
BACKSPACE 9 
BACKSPACE 9 
GO TO 32 

39 IF ( T DATUM ( K+l )-TDATUM (K+2I-TBEGIN) 40,40,38 

40 DO 41 KJ = 1 ,21 
K1 * K + KJ - 1 

41 DATUMT(KJ,J) - TDATUM(Kl) 

42 CONTINUE 

IF (MOON) 43,45,43 

43 DO 44 K J=4 , 2 1 

44 DATUMT1 K J, MOON) = DA TUMT ( K J , MOON ) *ERTOAU 

45 ASSIGN 46 TO NS1 

PART 10. PUT AWAY NEEDEO DATA. TEST NAME, TIME OF BEGIN AND END. DO NOT 
WRITE TAPE 3 UNTIL TBEGIN PREDATES THE END OF THE FITTED 
INTERVAL. 50 REPEATS OLD DATA, 57 WRITES NEW DATA. THE NAMES 
ARE ERASED FROM TMAKE AS SOON AS THE DATA POST DATES TEND. WHEN 
ALL NAMES ARE GONE, RETURN TD INPUT 300 TO SEE IF ANOTHER 
EPHEMER I S IS TO BE CONSTRUCTED. 

46 00 65 K= 1 , 232 , 2 1 
DO 47 J*l, KOUNT 

IF ( COMP ARF(TDATUM(K),TMAKE(J) ) ) 47,48,47 

47 CONTINUE 
GO TO 65 

48 SWT * TBEGIN-TDATUM(K+l)-TDATUM(K+2) 

IF (SWT) 49,49,52 

49 IF(KTAGIJ)) 50,52,50 

50 WRITE TAPE 3 , ( DA TUMT ( K J , J ) , KJ-1,21) 

52 DO 53 KJ=1,21 
K 1 = K + K J 

53 DATUMT ( K J, J) - TDATUM ( K 1- 1 ) 

IF ( J— MOON ) 56,54,56 

54 DO 55 KJ = 4,21 

55 OATUMT ( K J , J ) = DAT UMT ( K J , J ) • ERTOAU 

56 IF (SWT) 57,57,58 

57 WRITE TAPE 3 , ( DA TUMT ( K J , J > , K J= 1 , 2 1 ) 

58 IF(TEND-DATUMT(2, J)-DATUMT(3, J) > 59,59,65 

59 TMAKE(J) = 0. 

DO 60 KK=1, KOUNT 

IF ( TMAKE ( KK ) ) 65,60,65 

60 CONTINUE 

WRITE OUTPUT TAPE 6, 61, F I LE , TBEG I N , TEND, KOUNT, ( TMADE ( KK ) , 

IKK* l , KOUNT ) 

61 FORMAT ( 28H0EPHEMERIS COMPLETED, FILE=F3.,6H, FROM F10.1,3H TO 
1 F10.1, 4H FOR 12, 18H BODIES AS FOLLOWS / 12(2X,A6)J 

FILE * FILE ♦ l. 

END FILE 3 
GO TO 2 

63 WRITE TAPE 3, FILE 
REWIND 3 

REWIND 9 
TAPE3 * 3. 

DO 64 J=3, 4000 

64 C(J> = 0. 

RETURN 
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65 CONTINUE 
GO TO 32 

66 PRINT 67, T8EGIN 

WRITE OUTPUT TAPE 6,67,TBEGIN 

67 FORMAT 1 3 3H TBEGIN PREDATES 2437202.5, IT IS FIO.l) 

68 CONTINUE 
REWIND 9 
END 

REM BSFILE! I , J ) BACKSPACES TAPE I UNTIL IT IS POSITIONED JUST 
REM BEHIND THE J TH EOF MARK. 

REM 

ENTRY BSFILE 

PZE 

PZE 

PZE 

BCD IBSFILE 


BSFILE 

SXO 

*-4,1 


SXD 

*-4,2 


SXD 

*-4,4 


XEC* SITES) 


TSX 

StRER) ,4 


LXD 

BSFILE-2,4 


CLA* 

1,4 


TSX 

$( IOS) ,4 


CLA* 

$ ( RDS ) 


STA 

BSF 


ANA 

AO 7000 


STA 

BTT 1 


STA 

BTT2 


LXD 

BSFILE-2,4 


CAL 

2,4 


ANA 

=0777777700000 


ERA 

=0007400000000 


TNZ 

ONEARG 


CLA* 

2,4 


TZE 

BACK 


PDX 

.1 


AXC 

* + 1,4 


XEC* 

S( ICO) 

BTT1 

BTTA 

• • 


TRA 

*♦1 

BSF 

BSFA 

* * 


XEC* 

$ ( RDS) 


XEC* 

SIdSR) 


AXC 

* + l ,4 


XEC* 

S( TCO) 

BTT2 

BTTA 

• • 


TRA 

CHECK 


T I X 

BSF, 1,1 


XEC* 

$ ( RDS ) 

BACK 

AXC 

*♦1,4 


XEC* 

SI TCO) 


AXC 

*♦1,4 


XEC* 

S(TRC) 


NOP 



AXC 

* + 1,4 


XEC* 

SITEF) 


NOP 



LXO 

BSFILE-4,1 


LXD 

BSFILE— 3,2 


LXD 

BSFILE-2,4 


TRA 

3,4 

CHECK 

TXL 

BACK, 1, 1 


LXD 

BSFILE-2,4 


CLA 

ERR+1 


STO 

0 


CLA* 

1,4 


LDQ* 

2,4 

ERR 

TSX 

8 , 4 


TXI 

BACK, 0, 14 


PZE 

BSF ILE-2»0, ERR 

ONEARG 

CLA 

BSFILE-2 


ADD 

=01000000 


STO 

BSFILE-2 


LXD 

CHECK, 1 


TRA 

BTT 1—2 

A07000 

OCT 

7000 


END 




REM SKF ILE ( l , J J SKIPS TAPE I OVER J EOF MARKS. 
REM 


ENTRY 

SKF ILE 


PZE 



PZE 



PZE 



SKF I LE SXD 

*-3,1 


SKD 

*-3,2 


SXO 

*-3,4 


TSX 

S(RER) ,4 

CHECK LAST READ 

TEFA 

**1 


TEFB 

*♦1 


LXD 

SKF ILE— 1 , 4 


CLA* 

1,4 

PICK UP THE TPE NUMBER 

TSX 

U IDS) ,4 

SET UP THE TAPE ADDRESSES 

LXD 

SKF ILE— 1*4 

LOAD IT AGAIN — MAN 

CAL 

2,4 

IS THERE A SECOND ARGUMENT 

ANA 

=0777777700000 


ERA 

=0007400000000 


TNZ 

ONEARG 

NO SECOND ARGUMENT 

GQGO CLA* 

2,4 

PICK UP THE SECOND ARGUMENT 

TZE 

BUMP+l 

DID SOME DUMMY WANT NO FILES 

LOOP SUB 

=01000000 


RDS XEC* 

SI RDS) 

READ THE TAPE 

TCOA 

• 


TCDB 

• 


TEFA 

BUMP 

DID WE HIT 

TEFB 

BUMP 

AN ENO OF FILE 

TRA 

RDS 

GO READ SOME MORE 

BUMP TNZ 

LOOP 


LXD 

SKFILE-3, 1 


LXD 

SKF I LE— 2 , 2 


LXD 

SKF ILE- 1 , 4 


NOP 



TRCA 

* + l 

TURN OFF TAPE CHECK 

TRCB 

*♦1 


TRA 

3,4 


ONEARG CLA 

SKFILE-1 


AOD 

=01000000 

SET UP XR4 FOR PROPER RETURN 

STO 

SKF I LE- 1 


PXD 

0, ,0 

SET UP FOR ONE FILE 

TRA 

RDS 


END 



COUNT 

1200 



REM INPUT ROUTINE USING ARITHMETIC STATEMENTS. CF NASA TN D-1092 
LBL INPUT, 6 

ENTRY INPUT 

REM THIS IS SUBROUTINE INPUT. ITS CALLING SEQUENCE 

REM CONTAINS THREE ARGUMENTS AN IDENTIFICATION 

REM CODE NUMBER, THE FIRST LOCATION RELATIVE TO WHICH 
REM ALL DATA IS TO BE LOADED, AND THE FIRST LOCATION 
REM OF A TABLE TO BE USED BY THE ROUTINE. 

REM 

REM 

REM INCLUDED IN THIS ASSEMBLY ARE SUBROUTINES 


REM 

1 INPUT 


REM 

2 CHRCTR 


REM 

3 CLEAR 


REM 

4 COMPAR 


REM 

5 ERROR 


REM 

6 LOOK 


REM 

7 NAME 


REM 

8 NUMBR 


REM 

9 STORE 


REM 

10 TABLE 


REM 

11 TEST 


REM 

12 ACCUM, FIX, 

FLT, BINARY 

REM 

13 PRINX 


REM 

14 REAO. 


REM 



INTAPE PZE 

0, ,7 

LEWIS INPUT TAPE NOT STD. 

OUTAPE PZE 

0, ,6 

FORTRAN STANDARD OUTPUT TAPE 

1 NDX PZE 


STORAGE FOR IRA 

PZE 


IRB 

PZE 


IRC 

BCI 

It INPUT 


INPUT SXD 

INOX, 1 

SAVE INDEX REGISTER A. 

SXD 

INUX+1,2 

SAVE INDEX REGISTER B. 

SXD 

INDX+2,4 

SAVE INDEX REGISTER C. 

NZT* 

1.4 

IF THE IDENTIFICATION NUMBER IS Z 

TRA 

4,4 

RETURN TO THE CALLING PROGRAM. 

CLA 

= 1B35 


ADD 

2,4 

2,4 IS THE BASE LOCATION. 

STA 

SET 


STA 

LOCI 


STA 

L0C4 


CLA 

TSXBS 

OPEN BACKSPACE GATE 

STO* 

$ l l ink ) 

CALL CHAIN WILL BACKSPACE 

LOCA CLA 

I, A 

1,4 IS THE IDENTIFICATION NUMBER. 

STA 

NREG1 


AXT 

36,1 

INITIALIZE 36 

STZ 

1*1,1 

LOCATIONS 

T I X 

*-1,1.1 

TO ZERO. 

STO 

I LOCI 

MAKE NON • ZERO. 

CLA 

3,4 

3,4 IS THE LOCATION OF THE TABLE. 

STA 

LOCFC 

PREPARE 

STA 

NREG1-1 


AOD 

= 1B35 

THE 

STA 

LOCFA 

ARGUMENT STORAGES 

STA 

LOCKL 


TSX 

CLEAR, 4 

CLEAR THE VAR REGION. 


00020 

00030 

00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

00160 

00170 

00180 

00190 

00200 

00210 

00220 

00230 

00240 

00250 

00260 

00270 

00280 


00290 

00300 

00310 

00320 

00330 

00340 

00350 

00360 

00370 

00380 

00390 

00400 

00410 

00420 

00430 

00440 

00450 

00460 

00470 

00480 

00490 

00500 

00510 

00520 

00530 

00540 

00550 

00560 

00570 



L0CA1 

CLA 

=0076100000000 INHIBIT READING UNTIL 

00580 


STO 

READ. 

ARRAY RECORD REFRESHED 

00590 


AXT 

A3, 2 

A3 FORCES RECORD TO BE FILLED 

00600 


SXD 

1*2 

IN CHRCTR 

00610 


REM 

LOOK AT THE FIRST CHARACTER ON THE FIRST CARD 

00620 


REM 

IN SEARCH OF A $ 

SIGN. 

00630 

LOCA A 

TSX 

CHRCTR, 4 


00640 


SUB 

=H0000S0 

CHECK FOR A $ SIGN 

00650 

LOCA. 

STO 

WORO 


00660 


TSX 

COMPARE 


00670 


OCT 

24251 7630000 

D, E, FILE FLAG, T 

00680 


LXA 

NREG1 , A 

ZERO IF SD HAS BEEN READ. 

00690 


TXL 

*+2,4,0 


00700 


TXI 

*+1,2, A 

BEFORE SD ADD A TO INDEX 2. 

00710 


TXH 

ERRU , 2,7 

JUNK 

00720 


TXH 

SGNOUT ,2,6 

$17 BEFORE $D. FILE FLAG. OFF 

00730 


TXL 

*+3,2,5 

$E BEFORE $D 

00740 


TSX 

REAO.+l , A 

CRASH READ GATE 

00750 


TRA 

L0CA1 

SHOULD NOW HAVE $D CARD 

00760 


TXH 

LOCAD, 2, A 

FIRST SD. 

00770 


TXH 

LOCAJ, 2,3 

$T AFTER SD. 

00780 


TXH 

LOCCK ,2,2 

S 17 AFTER SD. FILE FLAG 

00790 


TXH 

LOLBG ,2,1 

SE AFTER SD. 

00800 


REM 



00810 

LOCAC 

LXA 

READ. , A 

SD AFTER SD. TEST IF BUFFER 

00820 


TXL 

ERRU, A, 0 

OVERWRITTEN 

00830 


REM 

THIS Ii THE PROGRAM RETURN. 

00840 

RTN 

LXO 

INUX, l 

RESET INDEX A. 

00850 


LXO 

INUX+1,2 

RESET INDEX B. 

00860 


LXD 

INUX+2 , A 

RESET INDEX C. 

00870 


TRA 

A, A 

RETURN TO CALLING PROGRAM. 

00880 


REM 

HUNT FOR THE = SIGN OF THE $ DATA CARD- 

00890 

LOCAD 

CLA 

=0076100000000 INHIBIT READING UNTIL 

00900 


STO 

REAU- 

SDATA FIELD SCANNED 

00910 


TSX 

CHRCTR, A 


00920 


TSX 

COMPAR, A 


00930 


BCI 

1, =00000 


00940 


TRA 

*+5,2,2 


00950 


TRA 

ERKO 

JUNK 

00960 


TRA 

LOCAD 

ALPHABETIC 

00970 


TRA 

ERKO 

NUMERIC 

00980 


SXD 

ALF , A 

= SIGN 

00990 


REM 

USE ALF MODE TO TEST ALL CHARACTERS. 

01000 


REM 



01010 


REM 

COMES HERE WHEN = 

= SIGN HAS BEEN FOUND. GET THE 

01020 


REM 

IDENTIFICATION NUMBER FROM THE CARD. 

01030 

LOCAF 

LXD 

I, A 


01040 


TXH 

*+2,4,43 


01050 


TXH 

LOCAG, A, A2 

CARD SCANNED OUT. 

01060 


TSX 

CHRCTR , A 


01070 


TSX 

COMPAR, A 


01080 


BCI 

1,,S+- 0 


01090 


TRA 

*+9,2,2 


01100 


TRA 

ERRM 

JUNK 

onio 


TRA 

ERRM 

ALPHABETIC 

01120 

LOCAE 

TSX 

BINARY, A 

FORM BIN WD IN VAR 

01130 


TRA 

LOCAF 

BLANK 

01140 


SXD 

ERSW, 2 

MINUS SET TO BY PASS. 

01150 


TRA 

LOCAF 

PLUS NO EFFECT. 

01L60 


STO 

SIGN 

DOLLARS 

01170 


REM 

COMES HERE TG CHECK THE REGION CODE AND THE 

01180 


REM 

VALUE APPEARING ON THE SDATA CARD. 

01190 

LOCAG 

CLA 

VAR 

COMMA 

01200 


TZE 

ERRU 

DATA SET NO. MISSING 

01210 


ALS 

18 


01220 


STD 

»• 

SAVE IDENT AT TABLE ( l ) , 

01230 

NREG1 

SUB 

»• 

PLACE FIRST ARG IN THIS ADDRESS. 

01240 


TNZ 

RTN 

0 IF CALL CODE = SDATA CDDE 

01250 


STZ 

ALF 

ALF = 0 MEANS NO ALF INFO. 

01260 


SXA 

NRfcGl , 0 

INDICATE SDATA IS READ. 

01270 


REM 

INST. BELOW ALSO 

EXECUTED AT READ., PLACED THERE BY CHRCTR 

01260 

TSXRD 

TSX 

READ. +1, A 

HERE SNEAK PAST READ. GATE 

01290 


SXD 

TESTJK , 0 


01300 


TRA 

LOCAN 


01310 


REM 



01320 


REM 

COMES HERE IF IT 

WAS A S TABLE CARD. 

01330 

LOCAJ 

TSX 

TABLE, A 


01340 


TRA 

L0CAN3 


01350 


REM 



01360 


REM 

COMES HERE IF AN 

ALPHABETIC CHARACTER WAS FOUND. 

01370 

LOCAK 

TSX 

NAME, A 


01380 


TNZ 

SET-1 

ZERO MEANS ON LEFT OF = SIGN. 

01390 


LXD 

JK1 , 1 

IF JK 1 DIDNOT INCREASE THEN 

01400 

TEST JK 

TXL 

ERRL , 1 , ** 

AN = SIGN WAS NOT USED. 

01410 


SXD 

TESTJK, 2 

SAVE JK 1 FOR NEXT TEST. 

01420 


CLA 

ILOC 

SAVE SIGN OF TABLE ENTRY. 

01430 


STO 

ILUC1 


01440 


TRA 

L0CAN2 


01450 


REM 



01460 


LXD 

JK, 2 

PREPARE TO ACCUMULATE THE NUM8ERS 

01470 

SET 

CLA 

* * , 2 

IN THE PSEUDO ACCUMULATOR. 

01480 


STO 

TEMP 


01490 


CLA 

1LUC 


01500 


TPL 

LOCAM 

MINUS MEANS FLOAT THE NUMBER. 

01510 


TSX 

FL T , A 


01520 


TRA 

LOCAM 


01530 


REM 



01540 


REM 

COMES HERE IF NUMERIC FIELD. 

01550 

LOCAL 

TSX 

NUMBER, A 


01560 


STO 

TEMP 


01570 



LOCAM 

TSX 

ACCUM»4 


ACCUMULATE RESULTS IN ACC. 

01580 


TSX 

CLEAR, A 



01590 


LXA 

WORD, A 



01600 


PXA 

0 , A 


+ WORD IN ACC FOR LOCAR 

0 1 6 1 0 


T XL 

LOCAR, A, 58 


NOT COMMA 

01620 


TXH 

LOCAR , A , 59 


NOT COMMA 

01630 


LXD 

JK1 ,2 


COMMA 

01640 


CLA 

ACC 



01650 


STZ 

ACC 


INITIALIZE 

01660 


LDQ 

I LOCI 


IS THIS VARIABLE FIXED POINT. 

0L670 


TOP 

LOCI 


NEGATIVE IS FIXED POINT. 

01680 


TSX 

FIX, A 



01690 

LOCI 

STO 

**»2 


STORE THE NUMBER RELATIVE TO BASE. 

01700 

LOCAN 

LXD 

JK 1,2 



01710 

L0CAN1 

TXI 

* + l ,2*1 


RAISE STORING INDEX BY ONE. 

01720 

LOC AN2 

sxo 

JK l » 2 


SAVE IT. 

01730 

LOCAN3 

LXD 

□ PER, 1 


ANY OPERATORS LEFT OVER. 

01740 


TKL 

*+J ,1,0 



01750 

ERRL 

TSX 

ERROR, A 



01760 


BCI 

1 • 0 ( L ) 



01770 


CLA 

ACC 


ANY DATA LEFT OVER. 

01780 


TNZ 

ERRL 



01790 


REM 




01800 


REM 

CALL THIS THE 

SWITCH HOUSE. 

01810 

LOCAO 

TSX 

CLEAR, A 



01820 

LOCAP 

TSX 

CHRCTR » A 



01830 

LOCAQ 

TSX 

COMPAR , A 



01840 


BCI 

1, . (0000 



01850 


TRA 

*+6,2,2 



01860 


TRA 

LOCAR 


SO, $T , OR OPERATORS. 

01870 


TRA 

LOCAK 


ALPHABETIC 

01880 


TRA 

LOCAL 


NUMERIC 

01890 


TRA 

LOCAT 


( SIGN 

01900 


TRA 

LOCAL 


DECIMAL 

01910 

IQCAR 

LXD 

OPER, I 


ANY OPERATORS LEFT OVER. 

Q 1920 


TXH 

ERRL, 1,0 


HIGH MEANS ALREADY HAS OPERATOR. 

01930 


SUB 

=HOQ0O$O 


SPLIT OFF $ FROM OTHERS 

01940 


TPL 

LOCA. 


IF + PROCESS t TYPE CH 

01950 


REM 

WHAT KIND OF OPERATOR IS THIS. 

01960 


TSX 

COMPAR, A 



01970 


BCI 

I , +— /* , o 



01980 


TXH 

ERRL, 2, 5 


REMOVE THE JUNK. 

0 L 990 


TXH 

LOCAN, 2, A 


COMMA 

02000 


SXD 

OPER, 2 


SAVE REST, WILL BRANCH IN SUB ACCU 

Q2010 


TRA 

LOCAP 


AFTER BOTH OPERANDS HAVE BEEN FOUN 

02020 


REM 




02030 


REM 

COMES HERE IF 

THE 

OCT OR ALF MODE. 

02040 

LOCAT 

TSX 

CHRCTR, A 



02050 


TSX 

COMPAR, A 



02060 


BCI 

1, JOAOOO 



02070 


TRA 

* + 5,2 



02080 


TRA 

ERRL 


JUNK 

02090 


TRA 

LQCAZ 


A CHARACTER 

U2100 


TRA 

LOCAU 


0 CHARACTER 

02110 


REM 




02120 


REM 

COMES HERE IF 

EMPTY PARENTHESIS WERE FOUNO. 

02130 


TSX 

CHRCTR, A 


) SIGN, GET NEXT CHARACTER. 

02140 


TOP 

* + 2 


MINUS FOR NEW CARO 

02150 


TSX 

TEST, A 


INSERT COMMA IF NEEDED. 

02160 


CLA 

ILOCI 



02 L 70 


STO 

ILOC 


PREPARE TO GET VALUE OF 

02180 


LXO 

JK 1 » 2 


CURRENT LEFT SIOE. 

02190 


TRA 

SET 



02200 


REM 

COMES HERE IF 

OCTAL M006 • 

02210 

LOCAU 

TSX 

CHRCTR, A 



02220 


SUB 

= HOOOQO > 



02230 


TNZ 

LOcAU 



02240 


TRA 

LOLAW 


) SIGN 

02250 

LOCAV 

LDQ 

VAR 



02260 


RQL 

3 


REPLACE TOP 3 BITS 

02270 


LGR 

3 


BY NEXT OCTAL CHARACTER 

02280 


ROL 

3 


PUT IN BOTTOM OF MQ 

02290 


STO 

VAR 



02300 


REM 

COMES HERE WHEN ) 

IS FOUND. 

02310 

locam 

TSX 

CHRCTR, A 



02320 


TOP 

*+2 


MINUS FOR NEW CARD 

02330 


TSX 

TEST, A 



02340 


LXA 

WORD, A 


CHARACTER TO IRC 

02350 


T XL 

LOCAV, A, 7 


OCTAL DIGITS 

02360 


T XL 

ERR J , A , 58 


ALPHABETIC, JUNK, B, 9. 

02370 


TXH 

ERRJ, A, 59 


SPLITS ( 

02380 

LOCAX 

LXO 

JK1 ,2 


COMMA 

02390 


CLA 

VAR 



02400 


TRA 

LOCI 



02410 


REM 

CONVERT THE NUMBER TO BINARY. 

02420 

LOCAV 

TSX 

BINARY, A 



02430 


REM 




02440 


REM 

COMES HERE IF 

ALF 

MODE. 

02450 

LQCAZ 

TSX 

CHRCTR, A 



02460 


TSX 

COMPAR, A 



02470 


BCI 

1,100000 



02480 


TRA 

*+5,2,2 



02490 


TRA 

ERRK 


JUNK 

02500 


TRA 

LOC AZ 


ALPHABETIC 

025 1 0 


TRA 

LQCAY 


NUMERIC 

02520 


REM 

COMES HERE WHEN ) 

IS FOUND 

02530 



LOCb A 

LXA 

VAR, I 

) SIGN 

02540 


TNX 

ERRK, 1,0 

ALF COUNT WAS ZERO. 

02550 


sxo 

ALF , I 


02560 


TSX 

CLEAR, 4 


02570 


TSX 

CHRCTR, 4 

PULL THROUGH CHARACTERS AND STORE 

02580 


SUB 

=017 

FILE FLAG, NEVER NEG. 

02590 


TZE 

ERRB 

COUNT WENT PAST E 0 JOB. 

02600 


TSX 

STORE, 4 

THEM ONE AT A TIME. 

02610 


T I X 

•”4 ,1,1 

GO BACK TILL NCHAR = 1 

02620 


LXD 

J, 1 


02630 


LXD 

MSH I FT , 4 


02640 


CAL 

BLANK 


02650 


LGR 

42,4 


02660 


ORS 

VAR + 1,1 

FILL IN PARTIAL WORD WITH BLANKS. 

02670 

L0CB8 

AXT 

1,4 

IRC TO 1 

02680 


LXD 

JKl , 2 


02690 


CLA 

J 

PREPARE TO STORE ALF WDS 

02700 


STD 

LOCBC 1 


02710 

locbc 

SXD 

JKl, 2 


02720 


CLA 

VAR+1,4 


02730 

LDC4 

STO 

**,2 


02740 


TXI 

* + 1,4,1 

J = J + 1 

02750 

LOCBC 1 

TXH 

LOCBD , 4, »* 


02760 


TXI 

LOCBC, 2,1 

JK= JK+ l 

02770 

LOCBD 

STZ 

ALF 


02780 


TSX 

CLEAR, 4 


02790 


TSX 

CHRCTR, 4 

LOOK AT NEXT CHARACTER. 

02800 


TOP 

*+2 

MINUS FOR NEW CARD 

02810 


TSX 

TEST, 4 

PUT IN COMMA IF NEEDED. 

02820 


SUB 

=HO0000, 


02830 


TZE 

LOCAN 

GO RAISE AND STORE JKl. 

02840 


REM 



02850 


REM THESE ARE ERROR 

CALLS 

02860 

ERRB 

TSX 

ERROR, 4 


02870 


BCI 

1 , 0 ( B ) 


02880 

ERRD 

TSX 

ERROR, 4 


02890 


BCI 

1 , 0 ( D ) 


02900 

ERR J 

TSX 

ERROR, 4 


02910 


BCI 

1,01 J) 


02920 

ERRK 

TSX 

ERROR, 4 


02930 


BCI 

1 , 0 ( K ) 


02940 

fcRRM 

TSX 

ERROR, 4 


02950 


BCI 

1 , 0 ( M ) 


02960 

ERRU 

TSX 

ERROR, 4 


02970 


BCI 

1,0(U> 


02980 


REM 



02990 


REM 

(E COMES 

HERE AFTER (D 

03000 

LOCBG 

CLA 

=0076100000000 NOP 

03010 


STO* 

((LINK) 

CLOSE BACKSPACE GATE 

03020 


TRA 

LOCAC 

RETURN 

03030 


REM PURPOSE OF (END 

CARD IS TO PROTECT FOR I EGN DATA FROM 

03040 


REM BACKSPACE WHEN 

CHAIN IS CALLED. 

03050 


REM 



03060 


REM 

END OF THE 

MAIN SEGMENT 

03070 


REM THIS A ROUTINE 

TO BACKSPACE THE INPUT TAPE WHEN A 

03080 


REM CALL CHAIN IS GOING TO SPILL THE BUFFER. 

03090 


REM THIS ROUTINE IS 

EXECUTED FROM CHAIN VIA THE ONE 

03100 


REM WORD SUBROUTINE 

(LINK) WHICH CONTAINS EITHER TSX OR NOP 

03110 

TSX6S 

TSX 

LOCBS, 4 

TO BE STORED AT (LINK 

03120 

LOCBS 

SXA 

* + 4,4 

SAVE INDEX 4 

03130 


CLA 

INTAPE 

INPUT TAPE NUMBER 



CALL 

(( 10S) 

SELECT INPUT TAPE 

03150 


XEC* 

(IbSR) 

BACKSPACE IT 

03160 


AXT 

* » , 4 

RESTORE INDEX 

03170 


TRA 

1,4 

RETURN TO THE CHAIN ROUTINE 

03180 


EJECT 



03190 


REM THIS IS SUBROUTINE CHRCTR. IT STORES SUCCESSIVE 

03200 


REM CHARACTERS FROM 

THE CARD AT LOCATION WORD, READS 

03210 


REM SUCCESSIVE CARDS INTO THE ARRAY RECORD, AND PRINTS 

03220 


REM (( 

TYPE CARDS. 

THE FIRST CHARACTER FROM A NEW CARD 

03230 


REM IS 

STORED IN WORD WITH A MINUS SIGN. 

03240 


REM 



03250 


REM 



03260 

CHRCTR 

SXD 

TEMP-10,2 


03270 


SXD 

TEMP-17,4 


03280 


LXD 

1,2 

CARD COL COUNT, SAW COUNT 

03290 


TXH 

*+2,2,83 

TOO EARLY TO READ. 

03300 


XEC 

READ. 

GATE MAY BE CLOSED 

03310 


LDQ 

Q 

HAS UNUSED CHARACTERS FROM BEFORE 

03320 


CLA 

SIGN 

ZERO OR ( GOES TO TAG 

03330 

LOCCA 

ALS 

6 

SHIFT LEFT 1 CHARACTER 

03340 


SLW 

TAG 

CLEARS OR PRELOADS TAG 

03350 

LOCCB 

LXD 

ALF, 4 

NONZERO MEANS ALF MODE. 

03360 

LOCCC 

TXH 

LOCCD, 2, 43 

SAW COUNT GIVES COL 81 = 43. 

03370 


TXH 

LOCCG ,2,42 

WAS COL 80 PROCESSED. 

03380 

LOCCD 

PXD 

0,0 

CLEAR ACCUMULATOR. 

03390 


LGL 

6 

SHIFT NEXT CHARACTER INTO ACC. 

03400 


T I X 

LOCCE, 2, 14 

COUNT DOWN BY 14 

03410 


LDQ 

RECORD+3, 2 

LOAD NEXT WORO 

03420 


TXI 

*+1,2,69 

JUMP BACK COUNTER. 

03430 

LOCCE 

TXH 

LOCCF , 4,0 

RETURN IF ALF MODE- 

03440 


PAX 

0,1 

MOVE CHR. INTO INDEX l 

03450 


TXH 

LOCCF, 1,48 

TRA MEANS GOOD CHARACTER. 

03460 


TXH 

LOCCC, 1,47 

TRA IF BLANK 

03470 


TXH 

LOCCF, 1,43 

TRA IF GOOD CHARACTER. 

03480 


TXL 

LOCCF, 1,42 

TRA IF GOOD CHARACTER. 

03490 


ZET 

TAG 

HERE ON ( 

03500 


TRA 

PRINT 

HERE ON (( GO PRINT 

03510 


TRA 

LOCCA 

( GOES TO TAG. 

03520 


REM 



03530 



LOCCF 

SXD 


1.2 

SAVE SAW COUNT 

03540 


STQ 


Q 

SAVE UNUSED CHARACTERS. 

03550 


ADD 


TAG 

ATTACH $ SIGN IF PRESENT. 

03560 


STO 


WORD 

SAVE THE CHARACTER AT WORD. 

03570 


LDQ 


SIGN 

SIGN OF MQ NEGATIVE IF NEW CARD. 

03580 


STZ 


SIGN 

CLEAR SIGN. 

03590 


STZ 


TAG 

CLEAR TAG OF ANY S 

03600 


LXD 


TEMP-17, 4 


03610 


LXD 


TEMP-10,2 


03620 


TRA 


1,4 

RETURN 

03630 


REM 

PRINT OUT THE 

i% CAROS. 

03640 

PRINT 

STQ 


U 


03650 


XEC* 


U TES) 

CHECK FOR QUIET BUFFERS. 

03660 


XEC 


READ. 

FETCH NEXT CARD. 

03670 


LDQ 


Q 


03680 


LGL 


6 

SPACE CONTROL SAFE IN ACC 

03690 


LDQ 


BLANK 


03700 


AXT 


4,4 

FILL END OF OUTPUT 

03710 


STQ 


OUT BUF + 19, 4 

BUFFER WITH BLANKS. 

03720 


T I X 


*— 1,4,1 


03730 


LGR 


6 

SPACE CONTROL BACK TO MQ. 

03740 


STQ 


OUT BUF 

STORE SPACE CONTROL. 

03750 


AXT 


14,4 


03760 


LDQ 


RELORD+2, 4 


03770 


STQ 


OUTBUF+ 15 , 4 


03780 


T I X 


*-2,4,1 


03790 


TSX 


PR I NX ,4 


03800 


TRA 


*+3 


03810 

LDCCG 

XEC 


READ. 

ALMOST ALWAYS A NOP. 

03820 


XEC* 


*< TES) 

WAIT FOR QUIET READ BUFFER. 

03830 


STZ 


TAG 

CLEAR THE $* CHARACTERS. 

03840 


AXT 


14,2 

FETCH CARD. 

0385Q 


LDQ 


INBUF+14 , 2 

14 WORDS 

03860 


STQ 


RECORD* 2 , 2 


03870 


T I X 


*-2,2,1 


03880 


CLA 


TSXRD 

OPEN READ. GATE 

03890 


STO 


READ. 


03900 

LOCCJ 

AXT 


84,2 

CARD COL 1 IS 84 

03910 


CCS 


= 0 

SET MINUS ZERO IN SIGN 

03920 


STO 


SIGN 


03930 


LGL 


12 

SAVE COLUMN 79 AND 80 



LDQ 


BLANK 

BLANK OUT COLUMN 01 TO 84 



LGR 


12 

MAY HAVE LOOK AHEAD 



STQ 


RECORD* 1 




LDQ 


RECORD-12 


03940 


TRA 


LOCCA 


03950 


REM 




03960 


REM 


COMES HERE ON END OF FILE FLAG 

03970 

LOCCK 

LXD 


TESTJK,4 


03980 


TXH 


R TN , 4, 0 

WAS DATA LOADED. YES RTN 

03990 

SGNOUT 

XEC* 


$ ( TES) 

WAIT FDR QUIET OUTPUT BUFFER 

04000 


AXT 


6,4 


04010 


LDQ 


OUT+6,4 


04020 


STQ 


0UTBUF*6, 4 


04030 


T I X 


*-2,4,1 


04040 


AXT 


13,4 


04050 


LDQ 


BLANK 


04060 


STQ 


OUTBUF+ 19 , 4 


04070 


TIX 


*- i , 4 , l 


04080 


TSX 


PR I NX ,4 


04090 


XEC* 


SITES) 

WAIT FOR QUIET BUFFER- 

04100 

LUCOUT 

CALL 


SEXIT 

THIS WAY OUT FOR KEEPS 

04110 


REM 




04120 

OUT 

0CI 


6 , 1END OF FILE INPUT TAPE JOB COMPLETE 

04130 


REM 




04150 


REM 

END 

OF THE SAP 

SUBROUTINE CHRCTR. 

04160 


EJECT 



04170 


REM 

THIS IS SUBROUTINE CLEAR. IT INITIALIZES 

04180 


REM 

NECESSARY PARAMETERS FOR SUBROUTINE STORE. 

04190 


REM 




04200 

CLEAR 

SXD 


J,0 

SET J TO 0. 

04210 


STZ 


VAR 

CLEAR VAR ( 1 ) • 

04220 


SXD 


MSHIFT.O 

RESET MSHIFT. 

04230 


TRA 


1,4 

RETURN TO CALLING PROGRAM 

04240 


REM 




04250 


REM 

END 

OF THE SAP 

SUBROUTINE CLEAR. 

04260 


REM 




04270 


REM 




04280 


REM 

THIS IS FUNCTION COMPAR. IT EXAMINES THE CURRENT 

04290 


REM 

CHARACTER AND TESTS IT AGAINST THE CHARACTERS 

04300 


REM 

FOUND IN THE ARGUMENT. ALPHABETIC AND NUMERIC 

04310 


REM 

SPLITS ARE MADE 

IF THE CHARACTER IS NOT FOUND 

04320 


REM 

IN 

THE ARGUMENT 

. THESE TESTS ARE COUNTED AND 

04330 


REM 

THE 

NUMBER LEFT 

IN INDEX 2 CORRESPONDS TO THE 

04340 


REM 

SUCCESSFUL TEST 

. IF NO TEST IS SUCCESSFUL 

04350 


REM 

THEN INDEX 2 CORRESPONDS TO THE TOTAL TESTS *1. 

04360 


REM 




04370 

COMPAR 

LDQ 


1,4 

USE FIRST ARGUMENT IN CALLING 

04380 


AXT 


1,2 


04390 

LOCDA 

PXD 


0,0 


04400 


LGL 


6 

PULL IN 1ST TEST CHARACTER. 

04410 


TZE 


LOCDD 

DONE IF ZERO. 

04420 


CAS 


WORD 

CHECK TEST WORD AGAINST CARD 

04430 


TXI 


LOCDA ,2,1 

CHARACTER. 

04440 


TRA 


LOGDC 

EQUAL. 

04450 

LQCDB 

TXI 


LOlDA, 2 , 1 

NOT EQUAL. GET NEXT TEST 

04460 

LOCDC 

CLA 


WORD 

CHARACTER. 

04470 


TRA 


2,4 

PROGRAM RETURN. 

04480 



LOCDO 

CLA 


2,4 

USE SECOND ARGUMENT IN THE CALLING 

04490 


POX 


Of 1 

SEQUENCE I DECREMENT) AS THE TEST 

04500 


TNX 


LOCDC. 1,1024 

FDR ALPHABETIC-NUMERIC SPLIT. 

04510 


sxo 


LOCDF ,1 

BECOMES INCREMENT 

04520 


LXA 


WORD, I 

CHARACTER TO IRA 

04530 


TXL 


LOCDC, 1,9 

NUMERIC 

04540 


TXH 


LOCDF, 1,57 

SPECIAL 0 ZONE, SX 

04550 


TXH 


LOCOE, 1,49 

ALPHABETIC 0 ZONE, NO / 

04560 


T I X 


*+2,1,32 

KNOCK OFF 11 ZONE EXCEPT - 

04570 


T I X 


*+1,1,16 

KNOCK OFF 12 ZONE EXCEPT ♦ 

04580 


REM 


+ ANO - SIGNS 

WILL BE 116)10, / WILL BE 117)10 

04590 


TXH 


LOCDF, 1,9 

SPECIAL 

04600 

LOCOE 

TXI 


LOCDF, 2,-1 

ADJUST IRB FOR ALPHABETIC 

04610 

LOCDF 

TXI 


L0C0C,2,*» 

ADJUST IRB FOR SPLIT 

04620 


REM 




04630 


REM 

ENO OF THE SAP SUBROUTINE COMPAR. 

04640 


EJECT 



04650 


REM 

THIS IS SUBROUTINE ERROR. IT IS CALLED IF AN 

04660 


REM 

ERROR HAS DETECTED ON ANY OF THE INPUT CARDS. 

04670 


REM 




04680 

ERROR 

SKA 


*♦2,4 

SAVE SOURCE 

04690 


XEC< 


SITES) 

WAIT FOR QUIET BUFFERS 

04700 


AXT 


**, 4 


04710 


CLA 


1,4 

GET PRINT ARGUMENT 

04720 


STO 


OUT BUF 


04730 


AXT 


1,1 


04740 


CAS 


R 


04750 


TRA 


*+3 

S THROUGH V 

04760 


TXI 


*+l,l,-l 

R 

04770 

MESSA 

PXD 


BLANK+4, 1 

A THROUGH N 

04780 


ANA 


= 7017 


04790 


ARS 


16 


04800 


ACL 


MEiSA 


04810 


STA 


102. 


04820 


AXT 


4,4 


04830 

102. 

LDQ 


* * , 4 


04840 


STQ 


OUTBUF+5,4 


04850 


T I X 


—2,4,1 


04860 


AXT 


14,4 


04870 


LDQ 


REC0R0+2, 4 


04880 


STQ 


OUTBUF+19,4 


04890 


T I X 


—2,4,1 


04900 


TSX 


PRINX.4 


04910 


XEC* 

SITES) 

WAIT FOR QUIET BUFFER 

04920 


AXT 


19,2 


04930 


CLA 


BLANK 


04940 


STO 


OUIBUF+19,2 


04950 


T I X 


*-1,2,1 


04960 


LOQ 


*H • 

PICK UP • 

04970 


LXO 


1,2 

SAW COUNT 

04980 


TXL 


*+2,2,71 

BACK UP IF OVER 71 

04990 


TXI 


*+3,2, -69 


05000 


RQL 


6 

ROTATE ACCORDING TO CHR PART. 

05010 


T I X 


*-1,2,14 

COUNT CHARACTER PART. 

05020 


STQ 


OUT 8UF+ 19 , 2 

STORE ACCOROING TO RESIDUAL 

05030 


TSX 


PR 1 NX , 4 

PRINT THE * 

05040 


XEC* 


SI TES) 

WAIT FOR THE • TO BE PRINTED 

05050 


LXD 


ERSW,4 

PICK UP ERROR SWITCH. 

05060 


TXL 


LOGOUT, 4,0 

NUN ZERO MEANS TRY NEXT SET 

05070 


AXT 


1208,4 

BYPASS MARK 

05080 


SXO 


BLANK, 4 

MARK BYPASSED CAROS 

05090 


LXA 


NRkGl ,4 

NONZERO IF THIS $DATA CARD. 

05100 


TXL 


*+2,4,0 


05110 


TSX 


READ. +1,4 

CRASH READ GATE IF SDATA CARD. 

05120 

L0CE8 

TSX 


CHKCTR , 4 

SKIP TO NEXT $DAT A AND TRY THAT SET. 

05130 


TQP 


LOCEB 


05140 


SUB 


= HQ 000$ D 


05150 


TNZ 


LOCEC 

TRA NOT A $DATA CARD 

05160 


STQ* 


NRbGl-1 

PUTS - SIGN IN TABLE(l) 

05170 


LXO 


BLANK+7 ,4 


05180 


SXO 


BLANK, 4 


05190 


LXO 


INDX+2,4 


05200 


TRA 


LOCA 


05210 

LOCEC 

AOO 


= 5 

TEST FOR END FILE FLAG 

05220 


TZE 


SGNOUT 

END FILE.. GET OFF 

05230 


TRA 


LOCEB 

OTHER 

05240 


REM 




05250 


REM 

ERROR MESSAGES. FIRST WORD ALSO USfcO AS A BLANK. 

05260 

BLANK 

BCI 


4, REDUNDANCY CHECK 

052 70 


BCI 


4, ILLEGAL CHARACTER 

05280 


BCI 


4, NO MANTISSA 

BEFORE E. 

05290 


BCI 


4, NO ENTRY IN 

TABLE 

05300 


BCI 


4 , ST YPE MISSING OR WRONG 

05310 


BCI 


4, E XPON. OUT OF RANGE 

05320 


REM 




05330 


REM 

END 

OF THE SAP SUBROUTINE ERROR. 

05340 


EJECT 



05350 


REM 

THIS IS SUBROUTINE 

LOOK. IT SEARCHES THE TABLE 

05360 


REM l 

FOR 

THE NAME STORED AT LOCATION VAR. IF FOUND, 

05370 


REM 

THE 

ACC IS NON-ZERO AT THE RETURN. 

05380 


REM 




05390 

LOOK 

SXD 


TEMP-12,4 

SAVE INDEX REGISTER C. 

05400 


CLA 


J 

SUBROUTINE. 

05410 


STD 


LOCFE 


05420 


AXT 


2,2 JK = 2 

IN INDEX 8 

05430 


AXT 


1,1 

J1 = 1 IN INDEX A 

05440 

LOCFA 

CAL 


** , 2 

CAL TABV(JK). 

05450 


T2E 


LOCFG 

NO ENTRY THIS VARIABLE 

05460 


STD 


LOCFO 

DECREMENT HAS NEXT 

05470 


ACL 


=0377777000000 


05480 


ANA 


=0377777000000 

ENTRY LOC. SAVE DECR 

05490 


SUB 


J 

ONLY. CHECK ENTRY LENGTH. 

05500 


TNZ 


LOCFD 

IF NOT THE SAME, LOOK AT NEXT ENTR 

05510 


PXD 


0,2 


05520 


PDX 


0,4 

JM = JK IN INDEX C. 

05530 


78 



L0CF6 

CLA 


VAR + 1 > I 

SEE IF VAR AND THIS 

05540 

LOCFC 

CAS 


* * * A 

ENTRY AGREE 

05550 


TRA 


*+2 

IF SO, CHECK REST OF NAME 

05560 


TXI 


*+2,4,1 

RAISE JM BY ONE. 

05570 

LOCFO 

TXI 


LOCFA-1,2,** 

IF NOT SO, GO TO NEXT ENTRY. 

055B0 


TXI 


*+1,1,1 

RAISE J1 BY ONE. 

05590 

LOCFE 

TXL 


LOCFB, 1,** 

FINISHED IF J1 IS GREATER THAN J. 

05600 


TSX 


CLEAR, 4 

CLEAR IF THE ENTRY AGREES. 

05610 

LOCFF 

CLA* 

LUCFA 


05620 


STD 


I LUC 

SAVE COMMON INDEX AT ILOC. 

05630 

LOCFG 

LXD 


TEMP-12,4 

PREPARE TO RETURN. 

05640 


TRA 


1,4 

RETURN TO THE CALLING PROGRAM. 

05650 


REM 




05660 


REH 




05670 


REM 

END 

OF THE SAP SUBROUTINE LOOK. 

05600 


EJECT 



05690 


REM 

THIS IS SUBROUTINE NAME, IT IS USED TO 

05700 


REM 

CORRELATE NAMES 1 

FROM INPUT CARDS WITH INTERNAL 

05710 


REM 

MEMORY LOCATIONS 

BY REFERRING TO THE TABLE. 

05720 


REM 




05730 


REM 




05740 

NAME 

SXD 


TEMP— 20 , 4 

SAVE INDEX C. 

05750 


REM 

GET 

THE REST OF 

THE VARIABLE NAME. STOP AT ANY 

05760 


REM 

NON 

ALPHANUMERIC 

CHARACTER. 

05770 

LOCGfl 

TSX 


STORE, 4 


05780 

LOCGC 

TSX 


CHRCTR,4 


05790 


TQP 


*+2 

MINUS FOR NEW CARD 

05800 


TSX 


TEST, 4 

COMMA MAY BE NEEDED. 

05810 


TNZ 


* + 3 

LOOK FOR ZERO. IF ZERO, MAKE IT 

05820 


ACL 


=HQOOOOO 

A LETTER 0 

05B30 


STQ 


WORD 


05840 

LOCGE 

TSX 


COMPAR, 4 


05850 


BCI 


1 , = I 0000 


05860 


TRA 


*+5,2,1 


05870 


TRA 


LOCGF 

JUNK OR OPERATORS 

05880 


TRA 


LOCGB 

NUMERIC OR ALPHABETIC 

05890 


TRA 


LOLGG 

t SIGN 

05900 


STZ 


ILOCI 

= SIGN 

05910 


REM 

GO 

TO THE TABLE 

LOOKUP ROUTINE IF AN - SIGN 

05920 


REM 

OR AN UPERATOR 

HAS FOUNO. 

05930 

LOCGF 

TSX 


LOOK, 4 

FIND THE NAHE IN TABLE. 

05940 


TZE 


ERRT 

NAME HAS FOUND IN TABLE IF NON— ZfcR 

05950 


LXA 


ILUC.2 


05960 


TRA 


LOCGL 


05970 


REM 




05980 


REM 

GO 

TO THE TABLE 

VARIABLE LOOKUP ROUTINE IF A 

05990 


REM 

( SIGN HAS FOUND 

. 

06000 

LOCGG 

TSX 


LOOK, 4 


06010 


TNZ 


LOCGJ 


06020 

ERRT 

TSX 


ERROR, 4 


06030 


BCI 


1.01T) 


06040 


REM 

CONVERT THE INDEX TO BINARY. 

06050 

LOCGH 

TSX 


BINARY, 4 


06060 


REM 

GET 

THE NUMERICS 

FOR THE INDEX TO THE VARIABLE. 

06070 

LOCGJ 

TSX 


CHRCTR » 4 


06080 


TXL 


LUCGH ,1,9 

NUMERIC 

06090 


TXL 


EKKC , 1,27 

JUNK 

06100 


TXH 


ERRC , 1 , 29 

JUNK 

06110 


TSX 


CHRCTR, 4 

I SIGN. GET NEXT CHARACTER. 

06120 


TQP 


* + 2 

MINUS FOR NEW CARD 

06130 


TSX 


TEST, 4 

COMMA MAYBE NEEDED. 

06140 


TSX 


COMPAR, 4 


06150 


BCI 


1. *00000 


06160 


TRA 


* + 4,2,1 


06170 


TRA 


LOCGK 

OPERATORS 

06180 


TRA 


ERRL 

ALPHABETIC AND NUMERIC 

06190 


STZ 


I LUC l 

- SIGN 

06200 


REM 




06210 

LOCGK 

CLA 


VAR 

COMPUTE STORING INDEX. 

06220 


ACL 


ILUC 


062 30 


PAX 


0,2 

STORE AODRESS AT DECREMENT WITHOUT 

06240 


TXI 


*+1,2,— 1 


06250 

LOCGL 

SXD 


JK, 2 

ACCUMULATOR OVERFLOW. 

06260 


CLA 


I LOCI 


06270 


LXD 


TEMP-20,4 

RESTORE INDEX C. 

06280 


TRA 


1.4 

RETURN TO CALLING PROGRAM. 

06290 


REM 


CONSTANTS AND ERROR CALL. 

06300 

ERRC 

TSX 


ERROR, 4 


06310 


BCI 


1,010 


06320 


REM 




06330 


REM 

END 

OF THE SAP SUBROUTINE NAME. 

06340 


EJECT 



06350 


REM 

THIS IS SUBROUTINE NUMBER. IT IS USED TO 

06360 


REM 

ASSEMBLE NUMERIC 

DATA FROM CARDS. ALL VALUES ARE 

06370 


REM 

TREATED AS FLOATING POINT NUMBERS IN THIS ROUTINE. 

06380 


REM 




06390 

NUMBER 

SXO 


TEMP-23,4 

SAVE INDEX C. 

06400 


SXD 


KNT2 ,4 

INITIALIZE 

06410 


STZ 


KNT3 

THE SUBROUTINE 

06420 


STZ 


KNTl 

BRANCH PARAMETERS. 

06430 


STZ 


KNT4 


06440 


STZ 


TEMP 


06450 


TRA 


LOCHB 


06460 


REM 




06470 

LUCHA 

TSX 


CHRCTR, 4 


06480 


TQP 


*+2, 

MINUS MEANS FROM NEW CARD 

06490 


TSX 


TEST, 4 


06500 


LOCHB 

TSX 

COMPAR, 4 


06510 


aci 

1..EOGOO 


06520 


TRA 

*+6,2,2 


06530 


TRA 

LOCHK 

JUNK OR AN OPERATOR 

06540 


TRA 

ERRE 

ALPHABETIC 

06550 


TRA 

LOCHC 

NUMERIC 

06560 


TRA 

LOCHE 

E 

06570 


CLA 

KNT2 

DECIMAL POINT. 

06580 


TNZ 

*♦3 

ZERO MEANS THIS IS THE SECOND POIN 

06590 


TSX 

ERROR , A 


06600 


BCI 

1«0(N) 


06610 


STZ 

KNT2 


06620 


STZ 

NEXP 


06630 


TRA 

LOCHA 


06640 

LOCHC 

CLA 

NEXP 

COUNT THE NUMBER OF DIGITS BEHIND 

06650 


ADD 

-1B35 

THE. IF THERE IS ONE 

06660 


S.TQ 

NEXP. 


06670 

LOCHO 

LXA 

KNTIt 1 


06680 


TXH 

L0CHD2 , 1 i 

10 DO NOT ACCUMULATE PAST 10 

06690 

LOCHOl 

TSX 

BINARY, 4 

CONVERT THE DIGIT TO BINARY. 

06700 


TZE 

LOCHA 

DO NOT COUNT LEADING ZEROS. 

06710 

L0CH02 

TXI 

•♦1.1*1 

COUNT TOTAL NO. OF DIGITS 

06720 


SXA 

KNT 1,1 


06730 


TRA 

LOCHA 


06740 


REH 

COMES HERE WHEN THE EXPONENT FIELD IS 

06750 

LOCHE 

CLA 

KNT1 

ENCOUNTERED. 

06760 


TNZ 

LOCHH 

THERE MUST BE AT LEAST ONE DIGIT 

06770 


TSX 

ERROR, 4 

BEFORE THE E OF AN E FORMAT NUMBER 

06780 


BCI 

l,0(S> 


06790 

LOCHF 

CLA 

KNT3 

SEE IF EXPONENT DIGITS HAVE ARRIVE 

06600 


TRA 

*♦2 


06810 

LOCHG 

CLS 

KNT 3 

SEE IF EXPONENT DIGITS HAVE ARRIVE 

06820 


TNZ 

LOCHK— 2 

NON ZERO MEANS SIGN IS OPERATOR. 

06830 


STO 

TEMP 

STORE SIGN OF EXPONENT. 

06840 


CLA 

KNT4 


06850 


TNZ 

ERRF 

NONZERO MEANS MORE THAN 1 EXP SIGN 

06860 


SXD 

KNT 4 , 2 

MAKE NOZERO. 

06870 

LOCHH 

TSX 

CHRCTR* 4 


06680 


TOP 

*+2, 

MINUS MEANS FROM NEW CARD 

06890 


TSX 

TEST, 4 


06900 


TSX 

COMPAR, 4 


06910 


BCI 

l,+-.000 


06920 


TRA 

*+7,2,2 


06930 


TRA 

LOCHK- 2 

OTHERS 

06940 


TRA 

ERRF 

ALPHABETIC 

06950 


TRA 

LOCHJ 

NUMERIC 

06960 


TRA 

ERRF 

DECIMAL 

06970 


TRA 

LOCHG 

MINUS 

06980 


TRA 

LOCHF 

PLUS 

06990 


REH 

CONVERT THE EXPONENT TO BINARY. 

07000 

LOCHJ 

CLA 

TEMP 


07010 


ALS 

2 


07020 


ADD 

TEMP 


07030 


ALS 

I 


07040 


ACL 

WORD 


07050 


STO 

TEMP 


07060 


SXD 

KNT3.2 

RECORD FACT FOR SECOND SIGN. 

07070 


TRA 

LOCHH 


07080 


REM 

COMES HERE WHEN AN OPERATOR WAS FOUND. 

07090 


CLA 

KNT3 

TEST FOR THE PRESENCE OF EXPONENT. 

07100 


TZE 

ERRF 

ZERO MEANS NO EXPONENT CAME. 

07110 

LOCHK 

CLA 

KNT2 


07120 


TZE 

*+2 


07130 


STZ 

NEXP 


07140 


CLA 

KNT I 

SEE IF MORE THAN TEN NUMBERS HAVE 

07150 


SUB 

-10B35 

BEEN CONVERTED 

07160 


TPL 

*+2 

IF SO, USE THE DIFFERENCE IN THE 

07170 


PXD 

0,0 

COMPUTATION OF THE EXPONENT. 

07180 


SUB 

NEXP 


07190 


ADD 

TEMP 


07200 


STO 

NEXP 


07210 


REH 

MANTISSA IN VAR AND THE EXPONENT IS IN NEXP. 

07220 


CLA 

VAR 


07230 


TZE 

LOCHQ 

SHORT CUT IF ZERO. 

07240 


LOQ 

*0233000000000 CHARACTERISTIC FOR LOW BITS 

07250 


LGR 

8 

LOW 8 BITS TO MQ 

07260 


RQL 

8 


07270 


LRS 

0 

BRING SIGN 

07280 


STQ 

VAR 


07290 


ORA 

-0243000000000 CHARACTERISTIC FOR HIGH BITS 

07300 


FAD 

VAR 


07310 


FRN 



07320 


STO 

VAR 


07330 


CLA 

NEXP 

THE EXPONENT 

07340 


A XT 

1.2 


07350 


LOQ 

-1. 

PUT 1 IN MQ 

07360 

LOCHL 

LBT 


EXPONENT IN ACCU 

07370 


TRA 

LOCUM 

FOUND NO BIT. 

07380 


TXH 

ERKV ,2,6 

EXPONENT EXCEEDS 63 

07390 


STO 

VAR-2 


07400 


FMP 

TAB+1,2 

THIS FORMS 10 **N6XP 

07410 


XCA 


SAVE IN MQ 

07420 


CLA 

VAR- 2 


07430 

LOCHH 

ARS 

I 


07440 


TZE 

LOCHN 

10**NEXP FINISHED. 

07450 


TXI 

LOCHL, 2,1 


07460 

LOCHN 

THI 

LOCHO 

MULTIPLY IF PLUS. 

07470 


FHP 

VAR 


07480 


FRN 



07490 


TRA 

LOCHQ 


07500 

LOCHO 

STQ 

VAR-2 


07510 


CLA 

VAR 

DIVIDE IF NEXP IS MINUS. 

07520 


FDP 

VAR-2 


07530 


XCA 


ANSWER BACK TO THE ACCUM 

07540 

LOCHQ 

LXD 

TEMP-23,4 

RESTORE INDEX C. 

07550 


TRA 

1,4 

RETURN TO CALLING PROGRAM. 

07560 


REM THESE ARE THE 

ERROR CALLS FOR SUB NUMBR. 

07570 


ERRE 

TSX 

ERROR, 4 


07580 


BCI 

worn 


□ 7590 

ERRF 

TSX 

ERROR , 4 


07600 


BCI 

1,0<F) 


07610 

ERRV 

TSX 

ERROR , 4 


07620 


BCI 

1,0(V) 


07630 


REH 



07640 


REM 

THIS IS THE FLOATING PT . TABLE USED 

IN DBC 

07650 


DEC 

1E+32, 1E+16, 1E + 8. 1E+4.1E+2 

CONVERSION TABLE 

07660 

TAB 

DEC 

1E+1 


07670 


REM 



07680 


REM 

END OF THE SAP SUBROUTINE NUMBER. 


07690 


EJECT 


07700 



REM 

THIS IS SUBROUTINE 

: STORE. IT STORES CHARACTERS 

07710 


REM 

AT THE ARRAY VAR. 


07720 


REH 




07730 

STORE 

SXD 

TEMP-13,1 


SAVE INDEX A. 

07740 


SXO 

TEMP- 14, 2 


SAVE INDEX B. 

07750 


LXD 

J* 1 


PUT J INTO INDEX REGISTER A. 

07760 

LQCJA 

LXD 

MSHIFT.2 


LOAD INDEX B WITH MSHIFT. 

07770 


T IX 

LOCJB, 2,6 


ADVANCE MSHIFT. 

07780 


AXT 

36,2 


REFRESH MSHIFT. 

07790 


TXI 

*+1,1,1 


RAISE J BY ONE IF MSHIFT IS DVER 

07800 


STZ 

VAK , 1 


CLEAR NEXT CELL 

07810 

LOCJB 

CLA 

WORD 


LEAVE SIGN BEHIND 

07820 


LOQ 

= 0 



07830 


LGR 

42,2 


MQVE CHARACTER 

07840 


STQ 

TEMP- 7 


PLACES TO THE LEFT. 

07850 


CAL 

TEMP-7 



07860 


ORS 

VAR+1, 1 


STORE THE CHARACTER AT VAR. 

07870 


SXD 

MSHIFT.2 


SAVE MSHIFT. 

07880 


SXD 

J.l 


SAVE J. 

07890 


LXD 

TEMP-13, 1 


RESTORE INDEX A. 

07900 


LXD 

TEMP-14,2 


RESTORE INDEX B. 

07910 


TRA 

1,4 


RETURN TO CALLING PROGRAM. 

07920 


REM 




07930 


REM 

END OF THE SAP 

SUBROUTINE STORE. 

07940 


EJECT 



07950 


REM 

THIS IS SUBROUTINE TABLE. IT IS USED TO 

07960 


REM 

CONSTRUCT A TABLE 

OF NAMES TO BE USED ON CARDS 

07970 


REM 

ANO THEIR MEMORY LOCATIONS RELATIVE TO ARG 2 OF 

07980 


REM 

THE CALLING SEQUENCE. 

07990 


REM 




08000 


REM 




08010 

TABLE 

SXD 

TEMP-15,4 


SAVE INDEX C. 

08020 


STZ 

TEMP 



00030 

LOCKA 

TSX 

CHRCTR, 4 



08040 


TSX 

COMPAR , 4 



08050 


BCI 

I, ,00000 



08060 


TRA 

* + 5 ,2 » 2 



08070 


TRA 

LOCKD+1 


JUNK 

08080 


TRA 

LOCKA 


ALPHABETIC 

08090 


TRA 

LOCKO+i 


NUMERIC 

08100 

LOCKE 

STZ 

TEMP 


COMMA 

08110 


TRA 

LOCKD 



08120 


REM 

COMES HERE TO 

CONVERT THE ADDRESS TO OCTAL FOR 

08130 

LQCKC 

CLA 

TEMP 


THE TABLE. 

08140 


ALS 

2 



08150 


ADD 

TEMP 



08160 


ALS 

1 



08170 


ACL 

WORD 



08180 


STO 

TEMP 


ADDS TO MAGNITUDE 

08190 


REM 




08200 


REM 

COMES HERE TO 

GET 

NUMERICS. 

08210 

LOCKD 

TSX 

CHRCTR, 4 



08220 


TSX 

COMPAR, 4 



08230 


BCI 

I , . =/000 



08240 


TRA 

*♦7,2,2 



08250 


TRA 

ERRA 


JUNK 

08260 


TRA 

ERKA 


ALPHABETIC 

08270 


TRA 

LOCKC 


NUMERIC 

08280 


TRA 

LOCKT 


/ CHARACTER 

08290 


TRA 

LOCKF 


- SIGN 

08300 


REM 




08310 


REM 

COMES HERE IF 

A DECIMAL PT WAS FOUND. 

08320 

LOCKE 

CAL 

TEMP 


DECIMAL POINT 

08330 


CHS 



MAKE SIGN MINUS 

08340 


STO 

TEMP 



08350 


TRA 

LOCKD 



08360 


REM 




08370 


REM 

COMES HERE IF 

AN 

- SIGN MAS FOUND. 

08380 

LOCKF 

TSX 

CLEAR, 4 



08390 

LOCKG 

TSX 

CHRCTR, 4 



08400 


TQP 

*+z. 


MINUS MEANS FROM NEW CARD 

08410 


TSX 

TESTT.4 



08420 


TZE 

LOCKH 



08430 


TSX 

COMPAR, 4 



08440 


BCI 

l,/,0000 



08450 


TRA 

*+5,2,1 



08460 


TRA 

ERRG 


JUNK 

08470 


TRA 

LOCKJ 


ALPHABETIC OR NUMERIC 

08480 


SXD 

B.2 


COMMA 

08490 


SXD 

a, 2 


SLASH 

08500 


TRA 

LOCKK 



08510 


REM 




06520 


REM 

COMES HERE TO 

STORE CHARACTER. 

08530 

LOCKH 

ACL 

-HOOOOOO 


REPLACE ZERO WITH LETTER 0 

08540 


STO 

WORD 



08550 

LOCKJ 

TSX 

STORE, 4 



08560 


TRA 

LOCKG 



08570 


REM 

COMES HERE AT 

END 

OF NAME. 

08580 



LOCKK 

TSX 

LOOK, A 


08590 


TNZ 

LOCKR 

GOES TQ LOCKR IF THERE IS AN ENTRY 

08600 


LXD 

J, 1 


08610 


TXI 

•+1.L.I 

ASSEMBLE KEY 

08620 


PXD 

0,1 

IRB HAS FIRST FREE LOC. 

08630 


ACL 

TEMP 


08640 

LOCKL SLW 

**,2 

STORE KEY INTO TABLE 

08650 


SXO 

*+1,1 

ADVANCE TO END 

08660 


TXI 

*+1,2,** 


08670 

LOCKM CAL 

VAR+1,1 

MOVE NAME, 0 TO TABLE 

08680 

LOCKN X£C 

LOCKL 

SLW IN TABLE 

08690 


TNX 

•+2,1,1 

TRANSFER WHEN DONE 

08700 


TXI 

LOCKM, 2,-1 

GO BACK TO FINISH 

08710 


ARS 

34 

KEEP ZONE OF 1ST VAR CH. 

08720 


TZE 

ERRG 

HAS NUMERIC, OR J=0 

08730 


REM 



08740 


REM 

REEXAMINE THE 

CUT OFF CHARACTER. 

08750 

LDCKP 

LXD 

B» 2 


08760 


TXH 

LGCKB ,2,1 

COMMA 

08770 

LOCKQ 

LXD 

TEMP-15,4 

/ CHARACTER 

08780 


TRA 

1,4 

RETURN. 

08790 


REM 



08800 


REM 

COMES HERE TO 

REPLACE KEY 

08810 

LOCKR 

ANA 

*0377777000000 J+l IN DECREMENT 

08820 


ACL 

TEMP 

LOCATION AND SIGN 

08830 

LOCKS 

XEC 

LOCKL 

SLW IN TABLE 

08840 


TRA 

LOCKP 


08850 

LOCKT 

CAL 

TEMP 

IS / LEGAL 

08860 


TZE 

LOCKti 

YES 

08Q70 


REM TRA 

ERRA 

NO, NUMERICS WAITING 

08880 


REM 



08890 


REM THESE ARE THE ERROR CALLS. 

08900 

ERRA 

TSX 

ERROR, 4 


08910 


BCI 

1 , 0 1 A ) 


08920 

ERRG 

TSX 

ERROR, 4 


08930 


BCI 

1,0(G) 


08940 


REM 



■ 08950 


REM END 

OF THE SAP SUBROUTINE TABLE 

08960 


EJECT 



08970 


REM THI 

S IS SUBROUTINE TEST. IT LOOKS AHEAD TO CLASSIFY 

08980 


REM A NEW CARD. ACOMMA WILL BE PUT INTO THE CURRENT 

08990 


REM CHARACTER POSITION ONLY IF EITHER ID THE NEXT 

09000 


REM CARD BEGINS WITH , 

A * SIGN FOLLOWED BY SOME OTHER 

09010 


REM CHARACTER OR C2I 

THE NEXT CARD BEGINS WITH AN 

09020 


REM ALPHABETIC AND AN 

= SIGN IS FOUND AND IT PRECEEDS 

09030 


REM 

ALL , $ AND . 

CHARACTERS ON THAT CARD. 

09040 


REM 



09050 


REM 



09060 

TEST 

SKD 

TEMP-12,4 

SAVE INDEX FOR RETURN. 

09070 


SUB 

-HOOOOSO 

TEST FOR A $ SIGN. 

09080 


TPL 

LOCLA 

POSITIVE MEANS $ SIGN. 

09090 


XEC 

READ. 

SAFE TO REFILL BUFFER 

09100 


T XL 

LOCLB, 1,16 

NUMBERS AND SPECIAL 

09110 


T IX 

*♦1,1,33 

FIX SO SLASH IS SPECIAL 



T I X 

• , i , 16 

MOD OUT ZONE 

09130 


TXH 

LOCLB, 1,9 

SPECIALS 

09140 


REM 


ALPHABETIC COME THRU. 

09150 

LOCLC 

AXT 

15,1 

SCAN THE CARO. 

09160 


LCQ 

RECORD+3, 1 


09170 


TXI 

*♦1,1,69 

FOR CHARACTER COUNT 

09180 


7 XL 

LOCLB, 1,70 

DONE IF WHOLE CARD SCANNED 

09190 

LOCLO 

PXO 

0,0 

OK TO -SEARCH 84 COLUNMS 

09200 


lgl 

6 


09210 


PAX 

0,2 

ZONE TO IRB 

09220 


ANA 

=017 

KEEP DIGIT 

09230 


SUB 

=013 

DIGIT PART OF ,*.= CHR. 

09240 


TZE 

LOCLJ 

CHECK ZONE 

09250 


T I X 

LOCLD ,1,14 

TRY NEXT CHARACTER 

09260 


TRA 

LOCLC+1 


09270 

LOCLJ 

TXH 

LOCLB, 2, 15 

, . S NEED NO COMMA 

09280 

LOCLA 

AXT 

B4,l 

84 IS CARD COL l 

09290 


SXD 

1,1 

RESET CHRCTR TO BEGIN CARD 

09300 


CLA 

RECORD-12 


09310 ’ 


STO 

Q 


09320 


CLA 

=H00000, 

SUBSTITUTE A COMMA. 

09330 


STO 

WORO 


09340 

LOCLB 

LXD 

TEMP-12,4 


09350 


CLA 

WORD 

IN AC FOR SR NAME, TABLE 

09360 


TRA 

1,4 

RETURN TO THE CALLING PROGRAM. 

09370 


REM 



09380 


REM 

FOR TABLE SUB STATEMENTS 


TESTT 

SXD 

TEMP-12,4 


09390 


AXT 

84,4 

IF NEXT CARD HAS VALID 

09400 


LDQ 

RECORD-12 

LEFT PART OF SUBSTATEMENT 

09410 

LOCN8 

PXD 

0,0 


09420 


LGL 

6 


09430 


PAX 

0,2 


09440 


TXH 

LOCLB, 2, 48 

0 ZONES EXCEPT BLANK 

094 50 


TXH 

LOCNC, 2, 47 

BLANK 

09460 


TXH 

LOCLB, 2, 27 

II ZONES AND ) 

09470 


TXH 

LOCLA, 2, 26 


09480 


TXH 

LOCLB, 2,11 

12 ZONES AND 8-4 

09490 


TXH 

LOCLA, 2, 10 

= 

09500 

LOCNC 

T I X 

LOCNB , 4 , 14 

NUMERICS AND BLANK 

09510 


LDQ 

RECORD+3, 4 


09520 


TNX 

LOCLB, 4,1 


09530 


TXI 

LOCNB ,4,70 


09540 


REM 



09550 


REM END 

OF THE SAP SUBROUTINE TEST. 

09560 


EJECT 



09570 


REM THE 

FOLLOWING FOUR 

SUBROUTINES ARE USED TO 

09580 


REM CONVERT DECIMAL DIGITS TO BINARY IN VAR, 

09590 


REM FIX 

FLOATING POINT 

NUMBERS, FLOAT FIXED POINT 

09600 


REM 

NUMBERS, AND 

FORM ARITHMETIC RESULTS IN THE 

09610 


REM 

PSEUDO ACCUMU 

LATOR [ACC) FOR EACH OPERATION 

09620 


REM 

ON A CARD. 


09630 


REM 



09640 


82 



BINARY CLA 

VAR 

ACCUMULATE A SERIES OF BASE 10 

09650 

ALS 

2 

DIGITS IN BINARY IN VAR. 

09660 

ADD 

VAR 


09670 

ALS 

L 


09680 

ACL 

WOKO 


09690 

STO 

VAR 


09700 

TRA 

1.4 


09710 

REM 



09720 

FLT CLA 

TEMP 

CONVERT TO FLOATING POINT THE 

09730 

LRS 

IB 

CONTENTS OF THE STORAGE CALLED 

09740 

ORA 

*0233000000000 TEMP. 

09750 

FAO 

*0233000000000 

09760 

STO 

TEMP 

LEAVE THE ANSWER IN TEMP. 

09770 

TRA 

1.4 


09780 

REM 



09790 

FIX UFA 

-0233000000000 CONV TO FIXED PT THE CONT 

09800 

LRS 

0 

OF THE ACCUMULATOR. 

09810 

ANA 

-0377777 


09820 

LLS 

0 


09830 

ALS 

18 

LEAVE THE FIXED POINT NUMBER IN 

09840 

TRA 

1.4 

THE ACCUMULATOR. 

09850 

REM 



09860 

ACCUN LXO 

0PER.2 

BRANCH FOR OPERATOR 

09870 

STZ 

OPER 

PREPARE FOR NEXT OPERATOR. 

09880 

CLA 

TEMP 


09890 

TRA 

*♦5,2 


09900 

TRA 

LOCMB 

• 

09910 

TRA 

LOCMA 

/ 

09920 

CHS 


MINUS 

09930 

FAD 

ACC 

PLUS 

09940 

ACCUN STO 

ACC 

NONE 

09950 

TRA 

1,4 


09960 

REM 



09970 

LQCNA CLA 

ACC 

01 V IDE. 

09980 

FDP 

TEMP 


09990 

STO 

ACC 


10000 

TRA 

1.4 


10010 

REM 



10020 

LOCMB LCQ 

ACC 

MULTIPLY. 

10030 

FMP 

TEMP 


10040 

TRA 

ACCUN 


10050 

REM 



10060 

REM 



10070 

REM END 

OF THE SAP SUBROUTINES ACCUH, FIX, FLOAT. 

10080 

EJECT 



10090 

REH SUBROUTINE PRINX 

DRAINS PRINT BUFFER TO LOGICAL TAPE 


REM 


GIVEN IN DECREMENT OF OUTAPE 


PRINX SKA 

PR I NY. 4 

BUFFERED WRITE ROUTINE 

10110 

CLA 

OUTAPE 

LOGICAL OUTPUT TAPE NUMBER 


CALL 

( IQS) 


10130 

AXC 

IOCD, 4 


10140 

XEC* 

A ( MRS) 

-WTDL 6 

10150 

XEC* 

A ( RCH ) 

-RCHL 0,4 

10160 

PKA 

0,4 

SAVE LOC OF 

10170 

STA* 

A ! MTC ) 

10 COMMAND FOR (HER) 

10180 

CLA 

TSXMR 

PRESET END ACTION 

10190 

STO* 

A(TES) 


10200 

PRINY AXT 

**,4 


10210 

TRA 

1,4 


10220 

T SXHR CALL 

(HER) 

EXECUTED FROM (TES) 

10230 

IOCD I ORT 

OUTBUF , , 19 


10240 

EJECT 



10250 

REM SUBROUTINE READ. 

FILLS READ BUFFER FROM LOGICAL TAPE 7. 

10260 

REAO. TSX 

READ. +1,4 

READ. GATE INITIALLY OPEN 

10270 

CLA 

-0076100000000 CLOSE READ GATE 

10280 

STD 

*-2 


10290 

SKA 

AXT, 4 


10300 

CLA 

INTAPE 

LOGICAL INPUT TAPE NO. 


CALL 

l IDS) 


10320 

AXC 

IOCD. ,4 


10330 

XEC* 

A(RDS) 

-RTDL5 

10340 

XEC* 

A I RCH) 

-RCHL 0,4 

10350 

CLA 

TSXTS 

SET UP BUFFER TEST 

10360 

STO* 

AITES) 


10370 

AXT ART 

• * , 4 


10380 

TRA 

1,4 


10390 

TSXTS TSX 

*♦1,4 

EXECUTED FROM (TES) 

10400 

SXA 

AXT, 4 

CLOSE OUT BUFFER 

10410 

CLA 

-0076100000000 SAY BUFFER IS QUIET 

10420 

STO* 

A(TES) 


10430 

AXT 

5,4 

PRESET REDUNDANCY 

10440 

SKA 

RTT, 4 

COUNT 

10450 

TSXTT AXC 

*+1,4 


10460 

XEC* 

AITCO) 

— TCOL 0,4 

10470 

AXC 

RTT, 4 


10460 

XEC* 

AITRC) 

-TRCL 0,4 

10490 

AKC 

X IT. ,4 


10500 

XEC* 

A(TEF) 

-TEFL 0,4 JOB COMPETE 

10510 

TRA 

AXT 

RETURN 

10520 

RTT AKT 

*• ,4 

INTERROGATE COUNT 

10530 

TIX 

SXA. ,4,1 

GIVE ANOTHER TRY 

10540 

AKT 

14,4 

CARD SURE BAD 

10550 

CLA 

INBUF+14,4 

SAVE IMAGE 

10560 

STO 

RECORO+2 , 4 


10570 

TIX 

*-2.4,1 


10580 

AXT 

64,4 

MAKE ERROR ROUTINE LOSE* 

10590 

SXD 

1.4 

IN INPUT BUFFER 

10600 

TSX 

ERROR, 4 


10610 

R BCI 

l.OIRl 


10620 

SXA. SXA 

RTT, 4 

SAVE COUNT 

10630 

AXC 

*♦2,4 


10640 

XEC* 

A(TEF) 

TURN OFF EOF IND 

10650 

XEC* 

A(BSR) 

BACKSPACE, 

10660 

AKC 

IOCO.,4 


10670 

XEC* 

A(RDS) 

REREAD 

10680 

XEC* 

A ( RCH ) 


10690 

TRA 

TSXTT 


10700 
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